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PREFACE 


This  report  is  one  of  a  group  documenting  National  Bureau  of  Standards  (NBS) 
research  and  analysis  efforts    in  developing  water  conservation  test 
methods,  analysis,  economics,  and  strategies  for  implementation  and  acceptance. 
This  work  is  sponsored  by  the  Department  of  Housing  and  Urban  Development/Office 
of  Policy  Development  and  Research,  Division  of  Energy  Building  Technology  and 
Standards,  under  HUD  Interagency  Agreement  H-48-78. 

The  report  was  prepared  by  Dr.  J.  A.  Swaf field,  Senior  Lecturer,  Drainage 
Research  Group,  Department  of  Building  Technology,  Brunei  University,  Uxbridge, 
UK.,  during  a  study  leave  period  as  a  guest  research  worker  at  NBS/Stevens 
Institute  of  Technology. 

Experimental  results  included  in  this  report  are  drawn  from  the  published  work 
of  the  Drainage  Research  Group  at  Brunei  University. 
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ABSTRACT 


The  flow  depth  and  velocity  changes  across  a  moving  solid  in  partially— filled 
pipe  flow  are  predicted  by  means  of  the  application  of  the  method  of 
characteristics  to  solve  the  unsteady  flow  equations. 

Simplified  force  models  are  presented  which,  when  used  in  conjunction  with 
empirical  relationships  linking  leakage  flow  past  the  solid  to  upstream  specific 
energy,  are  sufficient  to  provide  the  required  moving  solid  boundary  conditions 
that  allow  solid  velocity  prediction. 

A  wide  range  of  simulated  transport  conditions  are  presented  that  confirm  the 
applicability  of  this  technique  as  a  basis  for  the  future  evaluation  of  more 
complex  body  force  models. 

The  predicted  solid  velocity  during  drain  transport  is  shown  to  be  compatible 
with  laboratory  observations  of  the  influence  of  solid  dimensions  and  position 
in  inflow  profile  on  transport  characteristics. 

Key  words:     computer  based  model;  drainage;  solid  transport;  unsteady  flow 
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1.  INTRODUCTION 


The  prediction  of  solid  transport  in  partially  filled  drainage  pipe  flow 
requires  a  knowledge  of  the  forces  acting  on  the  solid  and  their  relative 
importance.     The  treatment  of  the  force  and  momentum  equations  for  the  body  and 
surrounding  flow  is  analytically  complex;  however,  a  numerical  model  capable  of 
simulating  the  unsteady  flow  equations  could  be  used  to  evaluate  models  of  the 
body  forces  for  comparison  to  laboratory  observations. 
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Previous  reports  [1,2]   have  developed  the  application  of  the  method  of 
characteristics  solution  to  the  unsteady  partially  filled  pipe  flow  equations. 
The  introduction  of  a  moving  solid  into  the  simulation  requires  the  knowledge 
of  leakage  flow  past  the  body  as  a  function  of  upstream  conditions,  with  a 
suitable  model  for  the  forces  acting  on  the  body.     The  passage  of  the  simulated 
solid  through  the  pipe  may  then  be  treated  as  a  moving  boundary  condition  by  a 
similar  technique  to  that  applicable  to  the  tracing  of  a  shock  wave  or  steep 
fronted  wave  motion. 

The  necessary  computational  techniques  to  allow  inclusion  of  the  solid  boundary 
conditions  are  presented  in  this  report  with  computer  simulations  involving 
three  force  models.     The  predicted  solid  motion  characteristics  are  also  com- 
pared to  laboratory  observed  solid  motion  in  drains  which  are  set  at  a  wide 
range  of  gradients. 

The  development  of  the  method  of  characteristics  solution  is  also  summarized, 
and  there  is  a  statement  of  the  necessary  boundary  conditions  at  the  drain 
inlet  and  exit  in  both  subcritical  and  supercritical  flow  regimes. 

It  is  stressed  that  the  objective  of  the  study  was  the  development  of  the 
necessary  computational  techniques  to  demonstrate  the  feasibility  of  the  method 
coupling  the  hydraulic  and  solid  motion.     The  development  illustrates  the  appli- 
cations of  the  basic  modeling  approach  and  the  simultaneous  solution  of  the 
governing  equations  of  hydraulic  interaction  with  that  of  the  motion  of  the 
solid.     A  basis  has  been  established  for  future  comparison  with  development  of 
other  body  force  models  representative  of  the  test  conditions  for  a  wide  range 
of  body  configurations  for  which  observed  laboratory  solid  transport  results 
are  obtained. 
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2.     DEVELOPMENT  OF  THE  NUMERICAL  SOLUTION 

2.1     SUMMARY  OF  THE  METHOD  OF  CHARACTERISTICS  SOLUTION  OF  THE  EQUATIONS  DEFINING 
UNSTEADY  FLOW  IN  PARTIALLY  FILLED  PIPE  FLOW 

Referring  to  figure  1,   it  has  been  shown  [1]    that  the  equations  of  motion  and 
continuity  may  be  expressed  as 

g  ill  +  g(S  -  S  )  +  V  II  +  11.    =     0  (1) 

8x  o         ax  at 
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VT-ii+  T  ill  +  A  ^    =  0. 
9x  8 1  3x 


(2) 


These  equations  may  be  combined  to  form  a  total  derivative  expression, 

dV  +  g  dh  +  g(g  _  s  )  =  0,  (3) 
dt       c  dt  o 

subject  to 

i2i  =  V  +  c  ,  (4) 
dt 

the  wave  speed  being  defined  by 


c  = 


gA 
T 


(5) 


where    V    =     local  average  flow  velocity 
c    =    local  wave  speed 
h    =     local  flow  depth 
Sq  =    pipe  slope 

2„2 

S    =     slope  of  the  local  energy  grade  line    =    ^  ^ —  ,  n  the  Manning 

m^/3 

coefficient  and  m  the  hydraulic  mean  depth 
T    =    water  surface  width 
A    =    water  depth  cross  sectional  area. 

Referring  to  figure  2,  if  the  variables  V  and  h  are  known  at  R  and  S  then  four 
equations  may  be  written  in  terms  of  the  unknowns  at  point  P,  by  means  of  a 
first  order  approximation, 

Vp  -        +  g_  (hp  -  h^)  +  g(S^  -  S^)ht  =    0  (6) 
^R 

xp  -  XR    =     (Vr  +  CR)At  (7) 

Vp  -  Vg  -  8_  (hp  -  hg)  +  g(S^  -  S  )At  =    0  (8) 

Xp   -  Xg     =      (Vg   -  Cg)At.  (9) 

It  is  stressed  that  these  equations  are  paired  and  that  equations  6  and  8  only 
apply  as  long  as  equation  7  and  9  are  satisfied.     This  introduces  the  charac- 
teristics lines  C"^  and  C~,  shown  in  figure  3. 

It  is  also  necessary  that  the  time  step  At  is  sufficiently  small  for  R  and  S 
(points  in  the  x-t  plane)  to  fall  within  +  Ax  of  point  P  as  shown  in  figure  2. 
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From  figure  3  it  will  be  seen,  provided  boundary  equations  governing  the 
conditions  at  the  extremities  of  the  system  are  known,   that  an  orderly  solution 
may  proceed  yielding  flow  depth  and  velocity  at  each  pipe  section  at  each  time 
increment. 

The  basis  for  the  system  boundary  conditions  is  set  out  below. 
2.2     PARTIALLY  FILLED  PIPE  FLOW  REGIMES 

Two  flow  regimes  may  be  identified  for  open  channels  or  partially  filled  pipes: 

(1)  Subcritical  flow,  Froude  N°    =  -^r  ^  ^ 

Here  the  local  wave  speed  exceeds  the  flow  average  velocity,  thus  waves 
may  be  propagated  both  upstream  and  downstream  in  the  flow,  i.e.,  c  >  V, 

(2)  Supercritical  flow,  Froude  N°  >  1 

Here  the  local  wave  speed  is  less  than  the  average  flow  velocity  at  that 
section  and  hence  waves  cannot  be  propagated  upstream,  i.e.,  V  >  c. 

The  flow  regime  applicable  to  any  partially  filled  pipe  flow  may  be  determined 
by  a  comparison  of  the  flow  normal  and  critical  depths. 

Under  steady  uniform  flow  conditions  the  force  balance  equation  for  an  element 
of  flow  is  normally  expressed  by  the  Chezy  equation 


m    =    hydraulic  mean  depth  A/P,  m  (where  A  is  the  pipeflow  cross  section 
area) 

Sq  =    sin  0  -duct  slope 
V    =    mean  velocity  m/ s 
C    =    Chezy  constant. 
P    =    pipe  wetted  perimeter. 

The  value  of  loss  coefficient  C  was  found  by  Manning  to  be  dependent  on 
hydraulic  mean  depth  and  duct  surface  roughness  n.     The  Manning  formula  is  the 
simplest  of  the  open  channel  equations: 


V 


C  /  m  Sq 


(10) 


V 


Q 


(11) 


where    Q  is  the  flow  rate  m-^/s 

A  is  the  flow  cross  sectional  area,  m' 
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The  value  of  the  Manning  coefficient,  n,  varies  with  pipe  or  channel  material; 
values  In  the  range  0,009  to  0,020  apply  to  for  materials  commonly  found  In 
building  drainage  systems. 


Equation  10  effectively  determines  the  flow  depth  under  steady,  uniform 
conditions  with  only  one  value  of  h  yielding  the  values  of  A  and  m  necessary 
to  satisfy  the  equation.     As  this  depth  Is  by  definition  constant  downstream, 
dh/dx  =0,  It  must  also  be  the  terminal  depth  corresponding  to  the  flow 
terminal  velocity  at  that  channel  slope. 

This  depth,  h^^,  Is  commonly  referred  to  as  the  normal  depth. 
The  specific  energy  of  the  flow  may  be  defined  as 

E     =    h  +  ^  (12) 

2g 

where    h    =     local  flow  depth,  m 

V    =     local  average  flow  velocity,  m/s 

From  equation  12  It  may  be  seen  that  the  flow  specific  energy  has  a  minimum 
value  below  which  the  given  flow  conditions  cannot  exist.     In  a  general,  non- 
rectangular  channel  this  value  may  be  determined: 

E     =  h+-^ 


2gA 
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^=0    =    1  ^  ,  (13) 

dh  g^3  dh 

From  figure  3 

dA  =  T  dh  (14) 
where      T  Is  the  surface  width  at  any  depth,  h. 

From  equations  (13)  and  (14)  the  minimum  value  of  E  will  occur  at  a  depth 
value,  he  that  satisfies  the  expression 

1  -  Q2T/gA3  =  0  (15) 
where      T  and  A  are  both  f (h) 

This  value  of  h  Is  referred  to  as  the  flow  critical  depth  h^. 

If  the  normal  flow  depth  h^^  exceeds  h^,  then  the  terminal  flow  would  be  termed 
subcrltlcal,  or  tranquil  flow.  If  h^  Is  less  than  h^  then  the  flow  Is  termed 
rapid  or  supercritical. 


6 


It  should  be  stressed  that  h^,  is  independent  of  pipe  slope  and  pipe  surface 
roughness,  while  the  normal  depth  is  dependent  on  both.     Thus  the  same  volume 
flow  rate  in  any  particular  pipe  may  be  rapid  or  tranquil  depending  on  pipe 
slope,  and  similarly,  the  same  flow  rate  in  a  series  of  constant  diameter  pipes 
will  be  tranquil  or  rapid  depending  on  roughness. 

Pipes  or  channels  in  which  rapid  flow  is  normal  are  termed  steep;   pipes  or 
channels  in  which  tranquil  flow  is  normal  are  termed  of  mild  slope. 

Figure  2  illustrates  the  importance  of  these  two  flow  regimes  on  the  solution 
of  equations  6  to  9. 

If  c  >  V,   then  the  conditions  at  P  are  determined  by  the  intersection  of  the 
C*"  and  C~  drawn  from  P  into  the  AC  and  BC  sections. 

If  c  <  V,   then  conditions  in  the  downstream  section  BC  cannot  affect  point  P. 
The  slope  of  the  C~  characteristic,  PS,  becomes  positive,  and  both  R  and  S 
lie  in  the  AC  section  as  shown. 

Both  subcritical  and  supercritical  flow  are  encountered  in  the  drainage 
applications  considered  and  the  equations  derived  below  cover  both  conditions. 

Referring  to  figure  2  for  subcritical  flow: 

Vr  *  OR) 


and 


-\  _ 

^C      ^R  _ 

^c 

-  ^A 

^C  "  ^A 

^c 

"  ^R  _ 

^C  ~  ^R  _ 

^c 

-  ^A 

^C  "  ^A 

he 

-  H  _ 

(Vr  +  cr)  - 

-  ^A 

At 
Ax 


as  xp    =    -KQ  and  xp  -  xr    =     (Vr  +  CR)At. 
The  solution  yields 

R       1  +  e  (V(,  -     +  cc  -  c^) 

.     c,(l  -  Vr  9)  +  C,V,  0  ^^^^ 
^  I  +         -  CpQ 

hR    =    he  -  (he  -  hA)(0(VR  +  cr))  (18) 

where      0  =At/Ax. 


Similarly, 
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i-e(v^-v, -cc+c,)  ^^^^ 

cs    =         11^^  ^""^  ^  (20) 

1+0   (cc  -  Cg) 

hg    =    he  +  0(Vs  -  cs)(hc  -  hg)  (21) 

For  the  supercritical  flow  regime,  the  equations  determining  R  in  figure  2 
remain  unchanged.     The  interpolation  equations  for  S'  in  figure  2  may  be 
determined  by  an  identical  technique  to  that  shown  above, 

Vj.  -  Vs«     =     (Vc  -  V^)  0(Vs'  -  cs') 

cc  -  CS'     =     (cc  -  ca)  0(Vs«  -  cs') 

and    he  ~  hs'     =     (h^  -  h^)  0(Vs'  -  cs')« 

The  solution  yields 

V^(l  -  c,0)  -  V,cp0 

Vet     =  ^  ^  ^  (22) 

l+0(Vc-V^  +  c^-  Cc) 

S  1  +  c^0  -  cc0 

From  equations  7  and  9  it  will  be  seen  that 

dt     =        1  (24) 
dx         V  +  c 

hence,  if  (V  +  c)  becomes  large,  then  At  becomes  small  for  a  constant  Ax.  This 
implies  that  the  progress  of  the  numerical  solution  could  become  prohibitively 
slow  for  supercritical  flow  conditions.     Fox  [3]    suggests  a  check  within  any 
program  that  will  terminate  the  solution  if  At  falls  below  a  specified  value; 
however,  this  comment  by  Fox  applies  to  the  existing  applications  of  the  method 
of  characteristics,  which  have  been  limited  to  large  civil  engineering  open 
channel  or  river  flooding  flow  problems.     It  is  likely  that  the  reduction  in 
At  will  not  be  a  significant  problem  with  the  relative  values  of  V  and  c 
encountered  in  drainage  sized  pipe  applications  due  to  Lheir  dependence  on 
pipe  geometry  and  flow  rate. 

The  determination  of  conditions  at  P  at  time  t  +  At  requires  the  following 
steps  (i-iv),  for  either  subcritical  or  supercritical  flow  conditions: 

(i)    All  conditions  known  at  time  t  for  nodal  points  A,  B,  C  etc. 
(figure  2). 
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(li)     Values  of  V,  h  and  c  at  Interpolation  points  R,  S,  or  S'  calculated 
from  equations  16  -  23. 

(iil)     Using  these  values  of  V,  h,  and  c,   the  conditions  at  P,  i.e. 

velocity  V  and  depth  h,   at  time  t  +  At,  are  calculated  by  means  of 
equations  6  and  8. 

(iv)     The  value  of  wave  speed  c  at  P  at  time  t  +  At  is  calculated  from 
equation  5.     The  value  of  flow  surface  width  and  cross  sectional 
area  are  calculated  from  flow  depth,  h,  and  the  channel  shape 
relationships . 

(v)     The  sequence  is  repeated  at  each  time  step. 

2.3     APPLICATION  OF  THE  METHOD  OF  CHARACTERISTICS  SOLUTION  TO  DRAINAGE  FLOW 

For  convenience,  the  application  of  the  solution  developed  above  to  drainage 
pipe  flow  may  be  considered  under  two  headings;   namely,  boundary  conditions  and 
characteristic  equation  solution  at  intermediate  pipe  sections.     Both  of  these 
headings  must  be  further  subdivided  depending  on  whether  the  flow  is  termed 
subcritical  or  supercritical. 

The  equations  6  to  9  may  be  restated  as 

Vp    =    X2  -  XI  h 

C+  (25) 

Xp    -    XR       =        (Vr    +  CR)At 

Vp    =    X4  +  X3  hp 

C-  (26) 

Xp  -  Xg    =     (VS  -  cs)At 
where      XI     =  g/cR 
X3    =  g/cs 

X2     =    Vj^  +  g  g(S^  -  S^)ht 

CR 

X4    =    Vg  -  g  ]^  -  g(S5  -  S^)At. 

It  will  be  noted  that  these  equations  apply  in  either  subcritical  or 
supercritical  flow,  the  interpolated  values  of  the  conditions  at  S  or  S '  being 
sufficient  to  define  the  flow  regime,  figure  2, 
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2.3.1     Initial  Steady  Conditions  Along  the  Pipe  Length 

As  will  be  seen  from  figure  3,  the  initial  conditions  along  the  whole  pipe  length 
at  time  t  =  0  must  be  known  in  order  for  the  solution  to  proceed.     It  is  there- 
fore necessary  to  calculate  the  steady  state  flow  velocity  and  depth  throughout 
the  pipe  length  initially  with  the  initial  wave  speed.     This  process  may  be 
carried  out  by  the  following  steps: 


(i)  Determine  the  steady  flow  normal  and  critical  depths  as  set  out 
previously.  This  determines  whether  the  flow  is  subcritical  or 
supercritical. 

(ii)     For  supercritical  flow,  the  normal  flow  depth  may  be  assumed  to 
apply  throughout  the  pipe  length.     As  the  velocity  exceeds  the 
wave  speed,  there  is  no  effect  propagated  upstream  from  the  pipe 
discharge  point.     This  implies  that  the  flow  leaves  the  pipe  at 
normal  depth  and,  for  a  known  flow  rate  and  pipe  dimension,  allows 
the  local  velocity  and  wave  speed  to  be  calculated  along  the  whole 
pipe  length. 

(iii)     For  subcritical  flow,  the  initial  water  surface  profile  is  more 
complicated  as  the  effect  of  the  pipe  discharge  is  propagated 
upstream.     In  subcritical  flow,  it  has  been  found  that  the  depth 
of  flow  is  at  its  critical  value  at  the  discharge  point,  with  the 
water  depth  then  rising  upstream  until  the  normal  steady  flow  depth 
is  achieved.     Calculation  of  this  depth  profile  is  presented  in 
[2]   and  summarized  below  in  terms  of  the  gradually  varied  flow 
profile  prediction  technique. 


Gradually  varied  flow  is  steady  non-uniform  flow  of  a  special  type.     The  flow 
parameters  are  assumed  to  change  slowly,  if  at  all,  in  the  flow  direction. 
The  basic  assumption  in  the  treatment  of  this  type  of  flow  is  that  the  local 
head  loss  at  any  section  is  given  by  the  Manning  expression  (11),  for  the  iden- 
tical local  flow  depth  and  rate  under  assumed  steady,  uniform  flow  conditions. 

Depth  profile  predictions  by  numerical  integration  are  based  on  this 
assumption, 

-;^{  ^+  (Zo  -  SoD  +  h  1      =  ^  (27) 

dL      2g  ^^2/3 

where       (Zq  -  SqL)  is  the  elevation  at  distance  L  along  the  channel,  measured 
in  the  downstream  direction;  Sq  is  sin  0,  channel  bed  slope. 


hence 


_  V  dV  +  s        dh     =     [  nQ  J  2 
g  dL        °      dL  Am2/3 


and  as,  Q  =  VA 
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dV  A  +  V  dA    =  0 

dL  dL 

and  as  =  T  from  equation  5,   it  follows  that 

dh 


dV    =     V  dA  =    -  VT  dh    =    _  QT  dh 
dL          A  dL  A  dL  A^  dL  ' 

and  substituting  in  equation  28  yields 

Q^Tdh  .   c  _  dh    ^    J     n  Q  i 
gA^  dL  dL  A 


(29) 


A 


dL    =    {     1  -  Q^T/g  a3  I 
S     -  (nQ/Am2/3)2 


hi 


1  -  Q^T/gA^ 


L     =    /   ^  

h      S    -  (nQ/Am2/3)2 
o  o 


dh 


(30) 


where      L  is  the  distance  between  two  known  depths  hg,  hi . 

The  initial  depth  at  each  section  Ax  apart  along  the  pipe  may  then  be 
calculated  from  the  profile  produced  by  the  integration  of  equation  30.  Flow 
velocity  is  then  calculated  based  on  a  constant  flow  rate  through  the  pipe,  and 
similarly,  wave  speed  is  determined  based  on  flow  depth  and  channel  geometry. 

Once  the  initial  flow  depth,  velocity  and  wave  speed  have  been  determined  the 
unsteady  flow  calculation  procedure  may  begin. 


2.3.2     Internal  or  Nodal  Points 


Simultaneous  solution  of  equations  25  and  26  at  all  points  Ax  apart  between 
X  =  Ax  and  x  =  (L  -  Ax)  will  yield  the  required  values  of  flow  depth  and 
velocity  at  the  end  of  each  time  step.     Wave  speed  may  then  also  be  determined 
from  equation  5.     This  process  applies  to  either  sub  or  supercritical  flow  con- 
ditions as  the  particular  regime  is  represented  in  the  interpolations  required 
to  fix  points  R,  S  or  S',  figure  2. 

2,3.3    Entry  Boundary  Conditions,  Supercritical  Flow 

In  this  case  the  inflow  profile  alone  determines  the  flow  depth  at  pipe  entry 
as  downstream  conditions,  that  would  have  been  represented  by  the  C~  character- 
istic in  subcritical  flow,  cannot  effect  the  flow  conditions  at  the  upstream 
boundary,  as  by  definition  the  flow  velocity  exceeds  the  wave  speed. 
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Hence  the  boundary  condition  is  obtained  from  the  flow  profile  Q  =  f(t)  solved 
with  the  normal  depth  expression 


Q    =     1  A  m2/3  s  1/2 
n  o 

where  A  and  m  are  both  f(h).     This  equation  may  be  rewritten  in  the  form 
2 

1  -  — (nQ)          =    0  (31 

a2  s 

o 


and  this  boundary  expression  may  be  solved  at  each  time  step  with  a  known  Q  by 
the  bisection  technique,  this  technique  is  described  later. 

2,3.4    Entry  Boundary  Conditions,  Subcritical  Flow 

For  subcritical  flow  the  downstream  conditions  do  contribute  to  the  entry  flow 
depth  and  velocity.     In  this  case  the  inflow  profile  Q  =  f(t)  is  solved  with 
the  C~  characteristic 


Q    =    f(t)     =  ViAi 
Vi     =    X4  +  X3  hi 

Q(t)     =    Ai   (X4  +  X3  hi)  (32) 

where  suffix  1  refers  to  the  entry  boundary  location  and  A  =  f(h)  dependent  on 
the  pipe  cross  sectional  geometry. 

In  the  form 

Q(t)  -  f(hi)(X4  +  X3  hi)     =    0  (33) 
this  boundary  equation  may  be  solved  by  the  bisection  technique. 
2.3.5     Exit  Boundary  Conditions,  Supercritical  Flow 

As  the  flow  velocity  exceeds  the  local  wave  speed  the  exit  boundary  condition 
may  be  determined  in  the  same  manner  as  the  upstream  nodal  points;  namely,  by 
the  simultaneous  solution  of  the  C^  and  C~  characteristics.     With  reference  to 
figure  2,  both  the  R  and  S'  points  lie  upstream  of  the  pipe  exit. 


The  exit  condition  may  be  included  in  the  nodal  point  calculations  for  the 
supercritical  flow  case,  no  separate  exit  subroutine  being  necessary  in  the 
program,  as  equations  25  and  26  are  sufficient. 
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2.3.6    Exit  Boundary  Conditions,  Subcritical  Flow 


At  pipe  exit  in  the  subcritical  flow  regime,   the  flow  depth  approaches  the 
critical  depth  value,  given  by  zero  value  of  equation  15: 

T     .       =  1 
 o   ^crit  ^» 

8 

crit 

where  A  and  T  are  f(h). 

This  condition  may  be  solved  with  the  C*"  characteristic 

V^+i     =    X2  -  XI  hfj+i 
where  N  =  N°  of  pipe  length  sections,  Ax. 
The  boundary  condition  becomes 

[(X2  -  XI  hN+i)  An+i1^  Tn+L   -1     =    0.  (34) 

8^N+1 

The  solution  may  again  be  achieved  by  use  of  the  bisection  method  with  the  use 
of  the  area  to  depth  relationship  for  the  channel. 

2.4     APPLICATION  OF  METHOD  OF  CHARACTERISTICS  TO  WASTE  SOLID  BOUNDARY, 
STATIONARY  CONDITION 

Considering  a  stationary  solid  deposited  at  some  point  along  the  waste  pipe, 
the  water  depth  and  velocity  upstream  of  the  solid  may  be  predicted  if  a  suit- 
able boundary  equation  may  be  written  linking  flow  past  the  solid  to  upstream 
conditions. 

Figure  4  illustrates  the  relationship  between  flow  past  a  stationary  solid  and 
the  specific  energy  upstream.     These  results  were  compiled  during  a  Brunei  Uni- 
versity Drainage  Research  Group  study  of  solid  transport  in  drainage  systems. 

The  flow  past  the  solid  (experimentally  determined)  may  be  expressed  as 

Q    =    K(h  +  SE„)2  (35) 

2g  o 

2 

where  SE  =  h  +  V  /2g,  flow  specific  energy  and  SE^  is  the  flow  specific  energy 
required  for  flow  initiation  past  the  solid. 

Equation  35  may  then  be  solved  with  the  CT*"  characteristic,  of  figure  5, 
Vx    =    X2  -  XI  hi 
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where  Q  =  Aj 


so  that 


At(X2  -  XI  hj)  =  K 


hj  +  —  (X2  -  XI  hj)^  - 
2g 


This  expression  results  in  a  quartic  in  terms  of  depth  upstream  of  the  solid, 
hx  (see  appendix  2).     This  quartic  must  be  solved  by  an  iterative  technique  as 
the  flow  area,  Aj^+i ,  is  a  function  of  hj.     The  Newton  Raphson  method  may  be 
employed  to  carry  out  the  necessary  iterative  solution.     Once  the  value  of  hj 
has  been  determined,   the  value  of        and  cj  may  be  determined. 

As  mentioned,  the  SEq  term  is  the  flow  specific  energy  required  to  initiate  flow 
past  the  stationary  solid.     If  the  value  of  flow  specific  energy  at  time  t  is 
less  than  that  of  SEq,  then  the  value  of  flow  velocity  at  the  solid  at  time 
t  +  At  is  set  equal  to  zero.     The  flow  depth  then  comes  directly  from  equation 


This  implies  that  the  flow  depth  upstream  of  the  solid  must  rise  to  SEq  prior 
to  the  initation  of  flow  past  the  solid.     This  solution  is  set  out  in  detail  in 
appendix  2,  and  it  has  been  shown  to  be  capable  of  simulating  depth  increase 
upstream  of  a  stationary  solid  [11. 

It  should  be  noted  that  the  analysis  above  applies  to  both  subcritical  and 
supercritical  flow  regimes,  as  only  the  C"^  characteristic  is  involved. 

2.5     APPLICATION  OF  METHOD  OF  CHARACTERISTICS  TO  WASTE  SOLID  BOUNDARY,  MOVING 
CONDITION 

Prior  to  application  of  the  moving  boundary  conditions  representing  the  solid 
motion,  it  is  necessary  to  determine  the  solid  motion  initiation  time.     This  may 
be  accomplished  by  monitoring  the  net  force  acting  on  the  solid. 

2.5.1     Model  of  Forces  Acting  on  the  Solid 

Figure  6  illustrates  the  forces  assumed  to  act  on  the  rectangular  section 
solid,  namely,  hydrostatic  pressure  forces,  Fpj^ ,  Fp2,  acting  on  the  trailing  and 
leading  edge  projected  areas,  body  weight  force,  frictional  force  due  to  wall 
to  body  contact,  Fg,  and  a  bouyancy  force,  Fg,  dependent  on  the  solid  position 
in  the  flow. 

The  net  force  acting  on  the  solid  may  be  expressed  as,  figure  6, 


25  as: 


hi     =  X2/X1. 


(36) 


Fpi  -  Fp2  +  mg  sine  -  Fg    =  FgoDY 


(37) 


where 


Fg    =     f(mg  cose  -  Fg) 
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and  f  is  the  wall  to  solid  friction  factor,   its  value  being  reduced  from  the 
static  friction  coefficient  to  the  sliding  friction  value  if  F30DY  becomes 
positive. 

If  FgoDY  ^  ^>  then  solid  motion  is  initiated  and  the  solid  velocity  at  the  end 
of  a  computational  time  step  At  is  given  by: 

where      AV    =    V|.+^^  ~  ^t         ™  ~  mass  of  body.  (39) 
(note  values  of  V^.  =  0  for  the  first  time  step) 

Hence.    AV    -    V^^^  =         ^t^W"  (*°> 

The  distance  traveled  by  the  solid  in  the  time  step  may  be  approximated  by 

=    X    +  1(V    ^  ^     +  V  )A  (41) 
B  o      2     t  +  At        t  t 

where    Xq  is  the  solid  position  at  time  t  and  Xg  the  final  position  at  the  end 
of  the  time  step. 

Equations  37  to  41  apply  to  subsequent  time  steps  with  the  mentioned 
modifications  to  the  value  of  wall  to  solid  friction  coefficient  in  equation 
37,     The  net  force  on  the  body  may  become  negative  in  subsequent  time  steps; 
however,  this  represents  solid  deceleration  and  no  modification  to  the  equations 
is  necessary  until  the  predicted  value  of  V^-  +        (equation  40)  becomes 
either  zero  or  negative  (the  condition  for  the  body  coming  to  rest). 

Figure  7  illustrates  the  forces  assumed  to  act  on  the  solid.     The  calculation 
of  the  buoyancy  force,  Fg,   (equation  37)  requires  the  solution  of  the  body  force 
system.     Values  of  Fj^,  F2  and  Fp  may  be  approximated  by  the  hydrostatic  equa- 
tion; however,  the  force  F^  is  not  readily  estimated  due  to  the  curvature  of  the 
solid  to  pipe  boundaries. 

Laboratory  observations  have  shown,  that  the  model  solids,  maternity  pads  [4]  , 
due  to  their  flexibility,  tend  to  take  up  the  shape  shown  in  figure  7,  Taking 
moments  about  the  leading  edge  point  A,  as  shown  in  figure  7,  allows  F^  to  be 
estimated. 

The  bouyancy  force  Fg  may  then  be  determined  as 

j^B    ^     ^^n  ~  ^d)  cosa  (42) 
where  a  is  the  slope  of  the  solid  surface  relative  to  the  pipe  wall. 
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Figure  8  illustrates  three  force  models  investigated  in  this  study: 


Model  1  -  buoyancy  forces  neglected,  the  value  of  Fg  =  0 

Model  2  -  downthrust  force  Fj)  =  0,     The  justification  for  this  model 

comes  from  test  observations.     In  the  boundary  equations  for 
the  solid  leading  edge,  i.e.  downstream,  the  leakage  flow  is 
assumed  to  take  up  its  normal  depth  at  once.     Laboratory  obser- 
vations indicate  a  downstream  transition  length  so  that  the 
flow  depth  immediately  ahead  of  the  solid  is  less  than  that 
predicted. 


Model  3  -  all  forces  shown  in  figure  7  included. 
2.5.2    Equations  Governing  Flow  Past  the  Moving  Solid 


It  may  be  assumed  that  the  relationship  between  specific  energy  and  flow  past 
the  solid  may  be  employed  in  the  moving  solid  case,  provided  that  the  absolute 
fluid  velocity  employed  in  equation  35  is  replaced  by  a  relative  water  to  solid 
velocity.     The  value  of  SEq,  the  specific  energy  term  to  initiate  flow  past  the 
body,  remains  unchanged. 

Equation  35  is  therefore  rewritten  as 
(V  -  V^)2  2 

Q     =    K(h  +   — ^  SE^)  (43) 

2g  o 

where      Q     is  the  leakage  past  the  moving  solid,  because 

.    Q    =    V^gA    =     (V^bs  -  Vb)A    =     (X2  -  Xlh  -  Vb)A  (44) 

where  A  is  the  flow  area  upstream  of  the  solid  and  V^^g  is  the  fluid  velocity 
expressed  by  the  C"*"  characteristic,  figure  5. 

Appendix  2  presents  the  full  solution  to  this  boundary  condition  in  a  general 
form,  applicable  to  both  the  initially  stationary  and  moving  solid. 

The  prediction  of  solid  velocity  allows  the  solid  path  to  be  drawn  in  the  x-t 
plane  as  shown  in  figure  9.     The  B"*"  lines  drawn  in  the  plane  are  the  equivalent 
to  the  fluid  characteristics;  the  gradient  of  B"*"  is  hence  given  by 

dx/dt    =    Vg.  (45) 

Figures  10  and  11  illustrate  the  necessary  techniques  to  allow  the  solution  to 
proceed  with  a  moving  solid.     A  slightly  different  solution  is  required, 
depending  on  whether  the  flow  is  subcritical  or  supercritical. 

Figure  10  presents  the  subcritical  case.    Assume  that  the  solid  was  at  point  C 
at  time  t  and  is  predicted  to  move  to  P'  by  time  t  +  At.     In  order  to  calculate 
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the  conditions  at  t  +  At,  It  is  necessary  to  set  up  the  interpolated  base  values 
at  R  and  S;  therefore,  for  conditions  immediately  upstream  and  downstream  of  the 
solid,  a  new  Ax  grid  A'P'   to  P'C  must  be  set  up. 

Using  the  position  P'  predicted  from  solid  velocity,   the  points  A'   and  C  are 
determined.     The  conditions  at  R  and  S  are  then  calculated  by  interpolation 
(equations  16-21). 

Point  R  may  be  used  to  yield  the  characteristic  RP '  that  may  be  solved  with 
the  solid  leakage  equation  43  to  yield  depth  and  flow  velocity  upstream  of  the 
solid. 

The  conditions  on  the  downstream  face  of  the  solid  are  calculated  by  a  similar 
technique  applied  to  C'B'.     The  C~  characteristic  is  solved  with  the  flow  rate 
at  the  solid  by  the  technique  utilized  at  pipe  entry  in  subcritical  flow 
(equations  32-33). 

Points  W,  X  and  Z,  figure  10,  may  be  dealt  with  by  the  nodal  point  equations, 
25  and  26,  as  the  necessary  interpolations  are  not  affected  by  the  presence  of 
the  solid  at  P ' . 

Conditions  at  P  and  Y,  however,   cannot  be  readily  calculated  due  to  the  B"*" 
characteristic.     Conditions  at  these  points,  however  are  required  as  base  con- 
ditions for  the  next  time  step.     As  the  Ax  and  At  values  are  small,  it  is  rea- 
sonable to  determine  conditions  at  P  and  Y  by  interpolation  between  X  and  P' 
and  P'  and  Z  respectively. 

Figure  11  illustrates  the  solution  for  the  supercritical  case.     It  remains 
necessary  to  approximate  conditions  at  P  and  Y  and  the  conditions  downstream  of 
the  solid  are  determined  by  the  use  of  equation  31   (the  normal  depth  calculation). 

2.6     SIMULATION  CASES  COVERED  BY  THE  PROPOSED  MODEL 

Two  types  of  solid  motion  must  be  covered  by  a  model  of  the  type  presented. 
The  motion  of  the  solid  subsequent  to  injection  into  the  flow  with  a  downstream 
velocity  must  be  dealt  with,  as  this  represents  the  introduction  of  solids  with 
finite  velocities  in  the  drain  from  water  closet  discharge.     Similarly,  the 
motion  of  a  deposited  solid  in  response  to  a  flow  must  be  considered,  as  this 
covers  the  partial  system  clearance  that  results  in  any  long  drainage  pipe 
exhibiting  a  series  of  depositions  along  its  length,  each  of  which  is  moved  on 
slightly  by  each  incoming  flow. 

Both  cases  may  be  dealt  with  by  the  proposed  model.     The  injection  case  is 
covered  by  assuming  that,  at  some  instant  of  time  and  position,  the  solid 
appears  in  the  flow,  moving  at  the  local  flow  velocity  that  was  calculated  at 
that  section  by  the  characteristics  solution  run  up  to  this  time  with  no  solid 
present . 
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The  deposited  solid  case,  at  any  downstream  location,  is  simply  dealt  with  hy 
the  technique  described  above  and  more  fully  in  [11.  In  this  case,  the  solid 
is  assumed  to  be  present  from  the  first  computing  time  step. 

Similarly,  the  model  must  be  capable  of  dealing  with  both  subcritical  and 
supercritical  flow  regimes.     The  techniques  for  these  cases  have  been  set  out; 
however,  care  must  be  taken  in  the  supercritical  case. 

Figure  12  illustrates  the  problem  encountered  in  supercritical  flow.     In  the 
case  of  a  stationary,  or  even  moving  solid,  the  depth  increase  upstream  of  the 
solid  may  be  sufficient  to  produce  a  subcritical  zone  in  an  otherwise  super- 
critical flow  condition.     In  order  for  the  upstream  supercritical  flow  to 
become  subcritical  behind  the  solid,  a  depth  change  wave,  or  series  of  depth 
change  waves  must  be  propagated  upstream.     If  the  depth  change  upstream  of 
the  solid  is  rapid,  as  would  happen  with  a  rapidly  increasing  inflow  profile, 
the  depth  change  waves  propagating  upstream  at  the  wave  speed  /g  A/T  combine 
to  form  a  steep  fronted  wave  that  moves  at  a  velocity  greater  than  /g  A/T,  and 
hence  the  solution  breaks  down. 

If  the  inflow  profile  is  not  sufficiently  rapid,  then  the  transition  from 
supercritical  to  subcritical  flow  may  be  accomodated  by  introducing  a  flow 
condition  check  into  the  calculations  at  the  stage  represented  by  figure  2.  If 
the  velocity  at  a  section  at  time  t  is  less  than  the  calculated  wave  speed  at 
that  section  in  an  initially  supercritical  flow  condition,  then  the  flow  locally 
is  assumed  to  have  become  subcritical  and  the  interpolated  values,   (figure  2) 
switch  from  R,  S'  to  R  and  S. 

The  speed  of  propagation  of  a  steep-fronted  wave  may  be  calculated  and  included 
in  the  model  in  a  manner  similar  to  that  described  for  the  solid  motion.  Fur- 
ther studies  will  include  this  facility.     The  test  cases  presented  in  this 
report  were  not  affected  by  this  effect,  due  to  the  inflow  profile  shape  and 
the  initial  position  of  the  solid  in  the  pipe. 
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3.     CALCULATION  TECHNIQUES  AND  PRESENTATION  OF  RESULTS 
3.1     DETERMINATION  OF  NORMAL  AND  CRITICAL  DEPTHS 

The  bisection  method  was  used  to  solve  the  equation  defining  both  critical  flow 
depth 

X    =     1  -  Q2T/gA3 
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and  normal  flow  depth 


Y    =  -  (n  Q/Am2/3)2 

that  feature  as  boundary  conditions. 

It  may  be  assumed  that  both  X  and  Y  have  zero  values  for  some  value  of  depth  h 
in  the  range  0  <  h  <  D  for  the  pipe  case. 

This  initial  interval,  0  <  h  <  D  is  bisected  and  h  =  D/ 2  used  to  evaluate  X,  Y. 
If  the  resulting  values  are  positive,  the  root  is  less  than  the  midpoint;  then, 
the  upper  limit  is  reset  equal  to  the  h  value  just  used  and  the  remaining 
interval  0  to  D/2  bisected.     The  process  is  repeated,  with  the  upper  limit 
replaced  until  the  value  of  X  and  Y  become  zero,  i.e.,  the  solution  desired 
for  the  critical  and  normal  flow  depths.     If  the  X  or  Y  value  had  been  negative 
then  the  root  would  be  greater  than  the  trial  h  value;   in  this  case  the  lower 
limit  would  be  intially  reset  to  the  trial  h  value,  D/s  (D/s  <  h  <  D) ,  and  the 
interval  bisected  to  a  new  value  of  3/4D,     The  process  is  repeated  until  a  root 
is  obtained. 

Due  to  the  need  to  include  the  area  depth  relationship,  this  solution  must  be 
undertaken  by  an  iterative  process.  The  time  taken  depends  on  the  complexity 
of  the  area-depth  function. 

3.2  NUMERICAL  INTEGRATION  FOR  SURFACE  PROFILES 
The  integration  of  the  position  vs  depth  profile 

h2 

L     =    /         1  -  QT^^/gA^  dh 
h^  S^  -  (nQAm2/3)2 

^1 

is  achieved  by  means  of  Simpson's  Rule.     Let  the  integral  X  =  /     F(h)  dh; 

^0 

then,  if  the  interval        -  hQ  is  divided  into  2  equal  increments,  the  value  of 
is  given  by 

X    =    ^  dh  [F(h^)  +  4F(h^  +  dh)  +  F(h„  +2  dh)]. 
2  o  o  o 

As  the  integration  moves  on,  the  length  traversed  may  be  accumulated  as  the 
added  interval  X  with  the  prior  L,  L  =  L  +  X,  at  the  completion  of  each 
integration. 

3.3  CHOICE  OF  TIME  STEP 

Referring  to  figure  4,  it  will  be  seen  that  the  time  step  chosen  must  be  such 
that  the  points  R  and  S  fall  within  +  Ax  of  point  P.     In  order  for  this  to 
occur  for  all  sections  along  the  pipe,  it  follows  that 
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This  expression  yields  the  smallest  possible  time  step,  as  the  maximum  values 
of  flow  velocity  and  wave  speed  at  any  pipe  section  have  been  used. 

In  order  to  ensure  that  the  computation  proceeds  as  quickly  as  possible,  the 
computer  program  presented  calculates  a  new  time  step  magnitude  at  each  time 
increment  so  that  the  time  step  increases  when  V  and  c  decrease,  but  decreases 
to  maintain  a  stable  solution  when  V  and  c  are  increasing  in  response  to  an 
inflow  surge. 

3.4  PRESENTATION  OF  RESULTS 

The  objectives  of  the  numerical  method  for  computation  of  the  transport 
mechanisms  of  solids  impartially  filled  pipe  flow  were  to:     (a)  identify  the 
potential  application  of  coupling  the  method  of  characteristics  solution  for 
the  hydraulic  phenomena  with  the  solution  of  the  equation  of  motion  for  the 
solid  based  upon  modeling  the  liquid/solid  interface  forces;  and  (b)  highlight 
any  limitations  inherent  to  the  technique.     The  simulated  pipe  flow/solids 
motion  numerical  data  were  developed  from  use  of  the  Fortran  computer  program 
TRANSCC,  run  on  the  NBS  CBT  Perkin  Elmer  732  computer. 

The  following  test  cases  were  investigated  and  are  reported  here 

(1)  Solid  types  —  two  cases  represented  by  figure  4,  270  mm  x  20  mm  x  70  mm 
and  270  mm  x  20  mm  x  35  ram  with  saturated  weights  of  250  g  and  125  g 

(2)  Pipe  gradients  —  1/40  to  1/300 

(3)  Inflow  profile  —  constant  profile  employed  with  peak  flow  of  4.2  il/s, 
overall  duration  9  seconds 

(4)  Pipe  roughness  coefficient  —  constant  at  0.015 

(5)  Solid  friction  factors  —  0.10  sliding  and  0.15  static 

(6)  Solid  position  in  flush,  3  s  and  4  s  from  flow  initiation  in  the  moving 
solid  insert  case,  at  2.7  m  from  entry 

(7)  Solid  position  —  4.2  m  from  entry  for  deposited  solid 

(8)  Force  models  —  three  models  as  illustrated  in  Figure  8, 

3.5  CHOICE  OF  CALCULATION  CONSTANTS 

In  order  to  undertake  the  calculation  procedure  described.  It  is  necessary  to 
have  values  of  leakage  constants  and  initial  specific  energy  for  the  solid. 
In  the  simulations  presented,  the  values  of  these  constants  are  drawn  from  the 
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results  illustrated  in  figure  4,  referring  to  transport  tests  at  Brunei  using 
large  deformable  solids. 

The  surface  to  pipe  wall  friction  factors  are  also  required.     No  values  are  at 
present  available  for  the  deformable  solids  tested;  however,  similar  tests  con- 
ducted in  the  Plumbing  Research  Laboratory  at  the  National  Bureau  of  Standards 
utilizing  impermeable  cylindrical  solids  suggest  values  for  friction  factors 
in  the  range  0.6-0.8  for  the  sliding  and  static  cases.     It  is  likely  that  the 
deformable  solids  will  have  lower  values,  due  to  the  presence  of  the  material 
saturating  water  that  will  tend  to  provide  a  degree  of  lubrication.     All  the 
simulations  presented  were  carried  out  with  values  of  0.1  sliding  friction 
factor  and  0.15  static  friction  coefficient. 
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4.     DISCUSSION  OF  UNSTEADY  FLOW  SIMULATION  RESULTS 

As  outlined  in  the  introduction,   the  objective  of  the  study  was  the  development 
of  a  technique,  based  on  the  method  of  characteristics,  that  would  allow  the 
transport  of  discrete  solids  in  partially  filled  pipe  flow  to  be  modeled 
mathematically. 

It  was  recognized  that  insufficient  data  on  the  frictional  characteristics  of 
model  solids  used  in  laboratory  test  programs  at  Brunei  and  NBS  were  available 
to  attempt  to  predict  actual  solid  transport  velocity  profiles.     However,  the 
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models  of  the  forces  acting  on  the  solid  were  based  on  current  research  at 
Brunei  and  NBS,  and  the  values  of  friction  factors  assumed  were,  as  already 
discussed,  based  on  the  best  available  Information. 

It  is  more  important  in  this  effort  to  demonstrate  the  ability  of  the  analysis 
techniques  to  yield  the  solid  velocity  and  the  associated  flow  depth  and  velo- 
cities during  the  transport  along  the  pipe.  Subsequently,  modifications  based 
upon  improved  force  models,  friction  factors  derived  from  tests  and  other 
empirical  adjustments,  can  be  introduced  to  fit  predicted  data  to  experimental 
results . 

From  the  examples  computed  for  the  conditions  of  3.4,  it  is  worth  noting  the 
form  of  the  relationships  of  solid  velocity  to  pipe  length  and  gradients. 
Figure  13  represents  typical  empirical  results  for  the  single  maternity  pad 
solid  tested  at  Brunei.     This  solid  has  dimensions  270  mm  x  20  mm  x  70  mm,  a 
saturated  weight  of  250  g,  and  a  specific  energy  to  leakage  flow  relationship 
in  a  100  mm  pipe  as  shown  in  figure  4.     This  solid,  with  a  half  width  version, 
was  used  as  a  basis  for  the  computer  simulations  presented. 

It  will  be  seen  from  figure  13  that  the  velocity  of  the  solid  over  the  major 
length  of  the  drain  is  governed  by  the  /L/G  term;  i.e.,  the  square  root  of 
distance  travelled  divided  by  pipe  gradient.     This  effect  will  be  studied  by 
plotting  the  computed  velocity  results  against  this  term.     Wherever  possible, 
simulations  are  presented  at  pipe  slopes  of  1/40  and  1/100  in  order  to  repre- 
sent both  supercritical  and  subcritical  flow  regimes.     In  addition,  these  gradi- 
ents represent  the  common  spread  of  slopes  employed  in  drainage  system  design. 

It  should  also  be  noted  that  the  water  depth  and  velocity  results  presented 
are  immediately  upstream  and  downstream  of  the  solid,  and  thus  refer  to  a 
location  that  moves  down  the  pipe  at  the  solid  velocity. 

4.1     INFLUENCE  OF  BODY  FORCE  MODEL  ON  PREDICTED  SOLID  VELOCITY 

Figure  14  illustrates  the  predicted  solid  velocity  for  the  three  force  models 
presented  in  figure  8.     It  was  found  from  the  values  of  force  predicted  by 
these  models  that,  for  the  range  of  model  sizes  simulated,  the  downthrust  on 
the  leading  surface  of  the  solid  materially  increased  the  surface  frictional 
force.     This  is  represented  in  figure  14  by  the  observation  that  at  both  1/40 
and  1/100  gradients,  the  model  3  results  yielded  the  lowest  solid  velocity.  As 
expected,  therefore,  the  exclusion  of  the  downsthrust  force,  Fq,  (figure  8) 
yielded  the  highest  solid  velocity.     The  omission  of  Fp  is  reasonable  on  the 
basis  of  laboratory  observations.     The  water  depth  tends  to  require  a  trans- 
ition length  downstream  of  the  solid  to  achieve  the  normal  depth  appropriate 
to  the  flow  past  the  solid.     Hence,  the  water  depth  predicted  by  the  analysis, 
the  effective  normal  depth,  will  be  an  overestimate  of  the  depth  of  this 
location. 

Figure  14  is  also  based  on  the  assumption  that  the  solid  appears  in  the  flow 
at  a  particular  x,t  coordinate  traveling  at  local  water  speed.     This  is  a 
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reasonable  assumption  that  leads  inevitably  to  a  rapid  deceleration  of  the 
solid  accompanied  by  a  buildup  of  water  depth  immediately  upstream  of  the 
solid.     The  adjustment  to  a  new  water-driven  velocity  is  rapid,  and  any  non- 
realistic  effects  due  to  the  insertion  model  are  likely  to  decay  rapidly. 

All  force  models  were,  used  in  the  collection  of  the  data  presented;  each  figure 
carries  a  note  as  to  the  applicable  model. 

4.2  INFLUENCE  OF  SOLID  DIMENSIONS  AND  POSITION  IN  THE  INFLOW  PROFILE  ON 
PREDICTED  SOLID  VELOCITIES 

Laboratory  experiments  and  data  analysis  presented  in  [4]  have  shown  that  the 
volume  available  behind  a  solid  in  the  flush  significantly  determines  its 
transport  properties. 

Referring  to  figure  15,  it  will  be  seen  that  introducing  the  solid  at  2.7  m 
from  entry  at  4  seconds  into  the  flush,  as  opposed  to  3  seconds,  does  support 
these  observations.     Although  the  predicted  differences  are  small,  they  are 
consistent. 

Similarly,  laboratory  tests  have  shown  that  base  area,  defined  in  the  case  of 
a  rectangular  model  solid  (figure  6),  as  thickness,  t,  times  width,  w,  is  a 
major  factor  in  determining  transport  performance.     Small  base  area  solids 
(e.g.,  tampons  or  sheets  of  toilet  tissue)  travel  badly  when  compared  to  larger 
solids.     This  is  due  to  the  increased  flow  past  the  solid  that  reduces  the 
depth  buildup  behind  the  solid,  and  from  equation  37,  leads  to  smaller  body 
forces.     Figure  16  clearly  demonstrates  this  effect  for  both  the  subcritical 
flow  at  1/100  and  the  supercritical  flow  at  1/40. 

4.3  COMPARISON  OF  SOLID  MOTION  FROM  BOTH  INITIAL  DEPOSITION  AND  INSERTION  AT 
WATER  VELOCITY  MODELS 

Figures  17  and  18,  and  19  and  20  compare  (for  slopes  of  1/40  and  1/100)  the  flow, 
solid  velocity,  and  depth  changes  for  a  solid  accelerated  from  rest  to  that  for 
a  solid  assumed  to  enter  the  flow  at  local  water  speed. 

A  number  of  general  observations  may  be  made  from  these  results  that  find 
confirmation  in  previous  observations  of  laboratory  tests: 

(1)  At  both  pipe  gradients,  the  maximum  depth  upstream  of  the  solid  occurs  in 
the  acceleration  from  rest  case. 

(2)  Perhaps  more  surprisingly,  the  maximum  depth  occurs  in  both  test  cases  in 
the  steeper  pipe.     This  is  a  direct  result  of  the  application  of 
equation  36. 

In  the  supercritical  flow  at  slope  1/40,  the  velocity  term  is  greater  than 
at  1/100,  and  consequently,  the  flow  depth  is  less.     Thus  the  "destruction" 
of  the  flow  momentum  by  the  partial  blockage  formed  by  the  solid  results 


25 


in  a  greater  depth  change  in  the  steeper  pipe.     Similarly,  this  enhanced 
depth  change  leads  to  a  higher  force  acting  on  the  solid  and  results  in 
the  earlier  motion  of  the  solid  in  the  1/40  pipe  case.     Figure  17  indi- 
cates solid  motion  at  1.4  s  for  1/40  and  4.2  s  at  1/100  pipe  gradients. 

(3)     Both  the  depth  upstream  of  the  solid  as  it  moves  along  the  pipe  and  the 
associated  solid  and  water  velocities  display  oscillations.     This  is 
entirely  attributable  to  the  choice  of  time  step,  and  the  fact  that  the 
simulation  equations  40,  41  are  only  first  approximations;  no  return 
iterative  technique  has  been  introduced.     Subsequent  time  steps  may  then 
underestimate  the  solid  velocity.     These  oscillations  are  reflected  in 
the  water  depth  upstream  of  the  solid.     This  link  between  upstream  depth 
and  solid  velocity  is  clearly  demonstrated  in  figure  18  by  the  sharp  dip 
motion  at  t  =  1.4  seconds.     All  simulations  involved  a  grid  length  of 
15/50  m,  time  step  dependent  on  wavespeed. 

4.4     COMPARISON  OF  PREDICTED  SOLID  VELOCITIES  TO  EMPIRICAL  RELATIONSHIPS 

Figure  13  indicated  that  solid  velocity  is  dependent  on  the  experimentally 
derived  group  /L/G, 

Vg  =   ci  -  C2  /l7g  , 

where  C^  and  C2  are  experimentally  measured. 

Figures  21  and  22  present  simulated  solid  velocity  results  plotted  against  the 
/L/G  term  for  two  force  models  and  pipe  slopes  from  1/40  to  1/300. 

The  _results  indicate  that  the  solid  velocity  in  each  case  is  linearly  dependent 
on  /L  over  the  major  portion  of  the  pipe  length.     The  dependence  on  G~l/2  is 
not  clearly  defined  by  the  results;  however,  it  is  clear  that  a  pipe  gradient 
terra  is  present  and  would  have  an  index  greater  than  0.5.     Alternatively,  and 
more  probably,  the  deviation  is  due  to  factors  not  yet  included  in  the  force 
models,  for  example,  the  surface  to  water  shear  force  arising  from  the  water 
flow  over  the  solid  has  not  been  included,  as  the  surface  shear  stress  is 
unknown.     However,  the  results  are  encouraging  in  that  the  general  form  of 
the  predicted  solid  velocity  curves  approximates  the  observed  laboratory 
results.     This  also  indicates  that  the  force  models  employed  were  reasonable, 
as  well  as  the  values  of  surface  to  wall  friction  coefficients  assumed. 

The  forces  acting  on  the  solid  have  not  been  presented  in  graphical  form,  as 
they  were  generally  found  to  remain  roughly  constant  during  any  simulation, 
although  exhibiting  the  oscillations  mentioned  previously.     Typical  values  for 
the  zero  downthrust  model,  figure  14,  at  1/40  were  -0.2  N  during  the  initial 
deceleration,   falling  to  -0.01  N  during  the  subsequent  motion.     At  1/100  the 
comparable  values  were  in  the  range  -0.2  N  initially,  falling  to  around  -0.02  N. 
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4.5     LIMITATIONS  TO  THE  SDIULATION 


As  fully  discussed  in  [2] ,  the  method  of  characteristics  solution  requires  an 
initial  flow  in  the  pipe  that  continues  beyond  the  termination  of  any  inflow 
profile.     This  effectively  means  that  a  simulated  solid  will  achieve  a  ter- 
minal velocity  in  the  pipe,  or  alternatively,  will  be  deposited  and  moved  on 
continuously  as  the  residual  flov/  acts  to  increase  upstream  depth.     This  effec 
is  not  readily  apparent  from  the  results  presented  due  to  both  the  relatively 
short  pipe  length,   15  m,  and  curtailed  run  duration  of  30  seconds.     For  the 
purpose  of  investigating  the  force  model  to  be  used,  this  limitation  is  not 
major  and,  indeed,  could  be  duplicated  in  any  parallel  experimental  work. 
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5.     CONCLUSIONS  AND  FURTHER  WORK 


The  objective  of  this  study  was  the  evaluation  of  the  coupling  of  the  method 
of  characteristics  with  the  equation  of  motion  of  a  solid  based  upon  an  assumed 
force  model  to  provide  a  numerical  analysis  base  for  modelling  the  transport 
of  discrete  solids  in  partially  filled  pipe  flow. 

The  results  presented  show  that  the  method  is  applicable,  and  that  the  motion  of 
the  solid  may  be  satisfactorily  tracked  through  an  x-t  grid  representing  the 
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pipe  by  introducing  a  solid  characteristic  whose  slope  in  the  plane  is  governed 
by  the  forces  acting  on  the  solid. 

A  wide  range  of  simulations  in  both  subcritical  and  supercritical  flow 
regimes  yielded  solid  velocity  results  that  were  compatible  with  laboratory 
observations . 

Further  work  is  required  to  establish  the  true  values  of  the  solid  sliding 
friction  coefficient,  with  further  study  of  the  force  models  to  be  used  as  the 
moving  solid  boundary  condition.     Similarly,  the  specific  energy  vs.  relative 
flow  rate  over  the  solid  requires  further  investigation  to  extend  the  range  of 
solid  geometry  available  at  present. 
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Figure  8.     Summary  of  alternative  force  models 
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Figure  21.     Predicted  solid  velocity  plotted  vs  /L/G,  pipe  slopes 
1/40  and  1/100 
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APPENDIX  I 
Program  TRANS CD 
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Program  TRANSCD 


This  appendix  presents  a  complete  print  out  of  this  program,  written  in  Fortran, 
a  flow  chart,  and  sample  input  data.     The  program  was  run  on  the  NBS  Center  for 
Building  Technology  Perkin  Elmer  732  computer. 

The  program  accepts  data  in  SI  units  with  the  exception  of  the  inflow  profile, 
which  is  read  in  as  litres/second  and  corrected  to  m^/s  within  the  program. 

The  program  is  written  in  terms  of  a  series  of  subroutines  that  deal  with 
specific  aspects  of  the  numerical  solution  presented  in  the  report.     In  order 
to  aid  the  understanding  of  the  program,  each  subroutine  is  discussed  individu- 
ally, with  the  calculation  methods  outlined  in  each  case.     The  titles  used 
refer  to  the  program  printout  included  in  this  appendix. 

Note  that  the  notation  employed  in  this  appendix  in  describing  the  various 
subroutines  is  compatible  with  that  used  in  the  main  report  and  not  the  program 
notation,  although  in  most  cases  the  variation  is  small. 

The  program  as  written  applies  only  to  simple  straight  pipes  with  choice  of 
constant  gradient,  diameter  and  roughness;  however,  these  subroutines  can 
obviously  be  utilized  in  the  modelling  of  more  complex  pipe  networks  in 
future  programs. 

Subroutine  TIMING 

This  subroutine  determines  the  maximum  values  of  wave  speed  and  velocity  at 
any  section  along  the  pipe  at  each  time  step.     This  ensures  that  the  next  time 
increment 

At  = 

(V  +  c) 

is  always  sufficiently  small  to  ensure  a  stable  solution.     This  effectively 
removes  the  need  for  a  time  step  factor,  but  this  is  retained  in  the  program 
and  normally  set  to  unity. 

Subroutine  LOSS 

This  subrutine  calculates  the  equivalent  steady  flow  friction  loss  over  each 
pipe  section  by  means  of  the  Manning  equation,  defined  as 

S    =    V  |V1  n2/ra^/3. 

Note  the  use  of  V  times  absolute  value  of  V  to  ensure  that  frictional  forces 
always  oppose  the  local  flow  direction. 
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Subroutine  DEPTH 


This  subroutine  is  only  called  at  the  initial  time  zero  to  calculate  the  normal 
and  critical  depths  based  on  the  initial  assumed  flow.     The  values  of  h^,  and  h^^ 
are  then  used  to  guide  the  solution  whenever  a  subcritical  or  supercritical 
flow  calculation  choice  arises. 

Subroutine  INFLOW 

This  subroutine  calculates  the  inflow  to  the  pipe  at  any  time  during  the 
simulation,  based  on  linear  interpolation  between  successive  pairs  of  0,  T 
coordinates . 

Subroutine  SHAPE 

This  subroutine  calculates  flow  surface  width,  T,  area  A,  and  wetted  perimeter 
depth. 

It  is  also  used,  at  the  initial  time  zero,  to  calculate  the  terms  necessary  to 
provide  the  subcritical  water  surface  profile. 

It  should  be  noted  that  this  subroutine  could  be  rewritten  for  any  other  pipe 
cross  sectional  shape  and  that  the  program  in  general  is  not  restricted  to 
circular  cross  section  pipes.     SHAPE  is  called  from  many  of  the  other 
subroutines. 

Subroutine  WAVSPD 

In  this  subroutine,  the  local  wavespeed  is  simply  based  on 


and  utilizes  SHAPE. 
Subroutine  PROFIL 

This  subroutine  is  only  fully  employed  if  the  initial  flow  is  shown  to  be 
subcritical  as  a  result  of  the  normal  and  critical  depths  calculated  via 
DEPTH. 

The  water  surface  profile  is  calculated  from  an  assumed  critical  depth  at 
pipe  exit  by  means  of  the  techniques  described  in  the  report.     The  calculation 
utilizes  SHAPE.     If  the  flow  is  supercritical,  then  PROFIL  sets  all  depths 
equal  to  the  normal  flow  depth. 


c 
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Subroutine  INTER 


This  subroutine  carries  out  the  necessary  interpolations  to  fix  conditions  at 
points  R  and  S  in  subcritical  flow  and  R  and  S'  in  supercritical  flow 
(figure  4). 

Subroutine  ENTRY 

This  subroutine  deals  with  the  flow  depth,  velocity  and  wave  speed  at  the 
upstream  boundary.     The  inflow  at  any  time  is  calculated  from  INFLOW  and  this 
flow  rate  is  then  employed,  in  the  subcritical  case,  in  conjunction  with  the 
appropriate  C~  characteristic  to  yield  the  required  parameter  values  referred 
to  in  figure  5. 

In  the  supercrititcal  flow  case,  identified  in  terms  of  the  initial  flow  normal 
to  critical  flow  depth  comparison,  the  inflow  curve  is  solved  in  conjunction 
with  the  normal  depth  relationship  as  the  downstream  conditions,  represented 
by  the  C~  characteristic  previously  employed,  can  no  longer  effect  conditions 
at  pipe  entry. 

Subroutine  NODAL 

This  subroutine  solves  the  CT*"  and  C~  characteristics  at  all  internal  pipe 
sections  between  x  =  Ax  and  x  =  L  -  Ax  (figure  4)  for  both  subcritical  and 
supercritical  flow. 

In  addition  for  supercritical  flow,  it  also  calculates  the  flow  exit  condition, 
which  is  based  solely  on  upstream  conditions. 

Subroutine  EXIT 

This  subroutine  is  only  utilized  for  the  subcritical  flow  case.     The  exit 
boundary  conditions  are  determined  by  solution  of  the  appropriate  C"*"  charac- 
teristic (figure  4)  with  the  critical  depth  expression 

g  a3 

and      Q     =    VA    =    A(X2  -  Xlh). 

As  this  is  effectively  a  loss  coefficient  concentrated  at  pipe  exit,  the 
boundary  equation  must  be  solved  in  the  form 

Q  IqI  T   _  1   =  0 

g  a3 

or      (X2  -  Xlh)   |(X2  -  Xlh)  |  =  1 

g  ^ 
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where  I    I   indicate  the  absolute  value  of  the  term  enclosed. 

The  solution  is  performed  by  means  of  the  bisection  technique  described  in  the 
main  text. 

Subroutine  ASSIGN 

This  subroutine  reassigns  the  values  calculated  at  the  end  of  a  time  step, 
points  P  in  figure  4  to  allow  the  calculation  procedures  to  move  on  one  time 
increment.     This  is  necessary  to  avoid  program  storage  space  problems. 

Subroutine  INSERT 

This  subroutine  ensures  that  the  solid  upstream  water  velocity,  depth,  and  wave 
speed  are  equated  to  those  calculated  at  the  insertion  time  and  pipe  section 
preset  in  the  input  data.     The  solid  velocity  is  set  equal  to  the  fluid  velo- 
city at  this  stage.     INSERT  is  only  called  once;  at  subsequent  time  steps, 
subroutine  SOLID  calculates  the  conditions  at  the  solid. 

Subroutine  SOLID 

This  subroutine  calculates  the  water  velocity  and  depth  upstream  of  the  solid 
by  means  of  the  specific  energy  vs.  relative  flow  rate  boundary  equation.  The 
flow  conditions  downstream  of  the  solid  are  determined  in  the  same  manner  as 
ENTRY,  dependent  on  flow  regime. 

Subroutine  FORCE 

This  subroutine  calculates  the  force  acting  on  the  solid.  Initially  it 

determines  the  inception  of  motion  by  monitoring  the  force  acting  on  the  solid. 

Various  force  models  may  be  investigated  by  modifying  this  routine. 

FLOW  CHART  PROGRAM  TRANS CD 

Set  up  initial  conditions. 

Time  =  0 

Read  pipe  data  -  Line  1  data  table 

Read  calculation  data  -  Line  2  data  table 

Read  inflow  profile  -  Lines  3-8  in  example  data 

Read  solid  position  -  Line  9  in  example  data 

Read  solid  dimensions  -  Line  10  in  example  data 

Read  specific  energy  to  flow  over  solid  coefficients  -  Line  11  in  example  data 
Read  force  model  indicators  and  friction  coefficients  -  Line  12  in  example  data 

Adjust  inflow  rate  from  lis  to  m-^/s  solid  dimensions  from  mm  to  m  solid 
saturated  weight  from  g  to  kg. 
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Calculation  of  initial  steady  conditions  at  time  zero. 


CALL  DEPTH  -  Calculation  of  normal  and  critical  depths  at  initial  flow. 


he  <  hn 

Subcritical  flow 
CALL  PROFIL  -  calculate 
water  surface  profile  and 
interpolate  for  section 
depths . 


^n  < 

Supercritical  flow 

CALL  PROFIL  -  set  all  depths 

to  normal  value. 


CALL  ASSIGN  -  set  up  base  arrays  H,  V,  C 

CALL  TIMINC  -  identify  maximum  V,  C  to  yield  stable  time  step. 
A"  =  PL/N 

At  =  Aoo/(V  +  c)n^x 

Output  -  Initial  conditions  and  pipe  description 

-  Solid  dimensions  and  model  indicators 

-  Depth,  velocity,  wavespeed  at  each  increment  if  N  <  10, 
at  1/10  pipe  length  points  if  N  >  10  and  a  multiple  of  10 

A  -  Update  time  and  commence  unsteady  flow  simulation. 

CALL  TIMINC  -  Identify  new  time  step  based  on  (V,  C)i^x 


Check  Time  vs.  TMAX 


Time  <  TMAX  Time  >  TMAX 

Go  to  B  Go  to  G 

B  -  Unsteady  flow  calculations. 


CALL  INFLOW  -  Calculate  pipe  inflow  rate  from  data  profile. 

CALL  INTER    -  Interpolate  to  obtain  base  conditions  HR,  HS  or  HS '  etc. 

Contains  choices  based  on  flow  regime  and  presence  of  the 
solid  beween  any  two  computing  points. 


CALL  LOSS      -  Calculation  of  equivalent  steady  flow  loss  terms,  SR,  SS. 

CALL  ENTRY    -  Solves  upstream  boundary  condition  for  either  sub  or  super- 
critical flow. 

CALL  NODAL     -  Solves  C+C~  characteristics  based  on  output  of  INTER  and  LOSS. 

Omits  the  two  nodes  bracketting  a  moving  solid  and  the  intial 
location  node  of  a  stationary  solid.     In  supercritical  flow, 
NODAL  also  provides  exit  conditions. 
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For  subcritical  flow  only,  CALL  EXIT  -  exit  depth  based  on  critical  depth  at 

pipe  discharges. 

Check  solid  model  being  simulated: 

Injected  solid,  time  <  injection  time,  Go  to  C 

Injected  solid,  time  _>  injection  time,  on  the  first  occasion,  Go  to  D 
Initially  stationary  solid.  Go  to  E 

Injected  solid,  time  >  injection  time,  following  first  occasion,  Go  to  E 

D  -  CALL  INSERT  -  assigns  water  velocity  at  solid  entry  station  and  time 
to  be  solid  velocity. 

Go  to  C 

E  -  CALL  SOLID  -  calculate  conditions  on  upstream  and  downstream  solid  faces. 
Downstream  velocity  given  by  same  technique  as  ENTRY  dependnet  on  flow 
regime. 

Go  to  C 

C  -  CALL  ASSIGN  -  sets  up  base  arrays  for  next  time  step. 
Output  of  simulation  results. 

Velocity,  depth,  flow  and  wavespread  at  solid  and  at  2  points  each  1/10 
pipe  length  upstream  and  downstream  of  the  moving  solid. 
Solid  velocity  and  body  forces. 
Solid  position. 

Check  presence  of  solid,  no  solid  Go  to  B,  solid  present  Go  to  F. 

F  -  CALL  FORCE  -  calculation  of  forces  acting  on  solid  over  next  time  step. 
Calculation  of  solid  position  and  velocity  at  end  of  next  time  step. 

Check  solid  still  in  pipe  at  end  of  next  time  step: 

Solid  position  >  pipe  length,  set  solid  indicators  to  no  solid  and  continue 
calculation  to  TMAX,  Go  to  B. 

G  -  Program  terminates. 

SAMPLE  INPUT  DATA,  PROGRAM  TRANS CD. 

Line  1. 


Pipe  diameter.  Manning  coeff.,  slope,  length. 
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Format  4F10.4. 

WW  -.1000  WW  0,0150  WW  0.0100  WW  15.0000 
Line  2. 

N°  calculation  steps,  Run  time,  Time  step  factor.     (Note,  N°  calculation 
sections  <  10  or  multiple  of  10  to  match  output  format,  time  step  factor, 
normally  set  =  1 . ) 

Format  13,  2F10.4 

V  50  WV  30.0000  WW  1.000 

Line  3. 

N°  pairs  of  coordinates  on  inflow-time  curve. 
Format  13 
W  5 

Line  4  to  8  (note  actual  number  depends  on  Line  3). 


Inflwo  QIN  or  time  TIN 
Format  2F10.4 
WW  0.2000  WW  0.0000 
WW  4.2000  WW  1.0000 
WW  1.2000  WW  4.5000 
WW  0.2000  WW  9.5000 
WW  0.2000  WV  50.0000 
Line  9. 

Solid  initial  position  in  pipe,  must 
location.     Anplie'?  to  both  initially 

Format  13 

V  10 


be  at  a  section  not  an  intermediate 
stationary  and  injected  models. 
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Line  10. 


Solid  dimensions,   length,  width,   thickness,  saturated  weight. 
Format  4F10.4 

VV  270.0000  VVV  70.0000  WW  20.0000  VV  250.0000. 
Line  11. 

Solid  minimum  specific  energy  to  allow  bypass  flow.   Specific  energy  vs.  Q 
coefficient. 

Format  2F10.4 

WW  0.0200  WW  0.6000 

Line  12. 

Buoyancy  force  indicator,  0.  force  omitted,   1.  force  included.  Downthrust 
force  indicator,  0.  force  omitted,  1.  force  included.     Solid  surface  shear 
stress,  set  0.     Sliding  friction  factor,  static  friction  coefficient. 

Format  5F10.4 

WW  1.0000  WW  0.0000  WW  0.0000  WW  0.1000  WW  0.1500 
Line  13. 

Solid  insertion  time,  zero  for  stationary  solid.     Solid  insertion  model 
indicator,  zero  for  stationary  solid,  >  0  for  moving  solid. 

Format  FlO.4  13 

WW  4.0000  W  4 

Units  -  SI  units  used  in  data  field  except  for  inflow  QIN  in  £.s,  converted  to 
m^/s  in  program.     Outut  in  S.I.,  except  flow  rate  in  Is.     Solid  dimensions 
and  saturated  weight  read  in  mm  and  g,  converted  to  m  and  kg  in  program. 
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tBATCH 

C  TRAMSC   ARf    A   SERIES  QF    PR0CRAH5  CFSIGNCO  TO  APPLY    TH£    METHOD  OF 

C  CHARACTERISTICS   SOLUriOH  TO  THE   EQUATIONS  DEflNlKC 

C  OKSTEAOY  FLOM    IN  fAi^TIALLY  FILLED  CHANNELS. 

C 

C 

01  HENS  ION  0P(61>,QIN( 30  I ,T INI  30 1 .VPf bl } ,HP 161 ) .CP ( 61) 
OIHENSION  V<61)*H(61I 

01  HENS  I  ON   Vf^(61),HRf6l),CR<61>«VS(61l«HSt6  1l,CS(6l  ).SS<61) 
DIMENSION   SRl6l>«X<^(6l)«XS(61l,XNtMI  ,C(61I 
COKMOH/CMl/LTG»UT»OX,D,SO*PL,SEO.XK,iM 
C0MK0N/CM2/hi THAX 
COMMON /Crt 3 /NPTS»OIN,riN 

COMMDN/CM'*/V,H,C,VR,HR  ,CR,XR,SP,  YS  ,HS,CS,XS.SS,XN 
COMH0N/CM5/HC  »HN 
COMMON /C M6/0P  » VP, HP tC  P 
C0MM0N/CM7/DL 

COMMON/CMfc/ INSOL,NSOL  l,NSOL 

COMHON/CM'5/SGLVELtHOS  ,VDS,CDS.HPLS,  VPOStCPDS 

COMMON/CMl  C/XPCS  XSnL,VPS0L,*SrL,O«»SUL.HuStVJS»CuS7HPUS,VPUS,CPUS 
COMMON/CMl 1/TL,Th,TT, SH 

COMMON/CMli'/Xf  AC2,  XE  AC,TAU,f  MC  V,f  STAT 

c 

C 
C 

c 

C 
C 

C  SETINITUL  CCNCIT  fQr<S 

TlME-0.0 

VPSOL '0.0 
C  READ  PIPE    DESCklPrns  DATA 

READ!  'itlOC  )r  ,  PKtSO 
C  0-DlA.    PM-HANNIMC   C1EFE.   SO-PlFt    SLOPE.    PL-PIPf  LfcNGlH 

100         FORHAK^FIC.',  ) 
C  READ  CALCLLAl  1  OH  0*1^ 

READ(^,103)^■»1^AX,TF4C 

IF(N.LE.IO)  IS2-I 

IF(N.CT.IO)  N2-N/10 
C  HZ    INCREflfhT    fiETH'-N   OUTPLT  hCTES.ALLOhS    nORc    THA^  lu 

C  CALCULATICh  NCCES    TT    BE  USED. 

C  N-NUMRFR    OF    PIPE    ScCriQNS   CChSICEPED,    T  MA  x-Oj  k  AT  1  0'-   OF  CalC. 

C  TFAC-TIhE   STEP  FA:n-<,l-lC. 

103         FORMAT! I3,2F1C.A) 

C  k£AO    INFLCV  LATA  P-^QFILf.    INfLCW  PRPrlLE    USti)    IS    dASfL.  jN 

C  A  LINEAR    INTE  PPCLATI  )N    3ETWEL^    Tr-FSE    DATA   POINTS.    NOIt  QftTA 

C  READ    IN    I^.  L/S   BUT  USED  FfOCEAH    I«  M«*J/S. 

RE  A0(  ^, 102 )h?  TS 
102  F0RMATU3) 

00   50  I-1,NPTS 

ftEAD(^«10A)01N(I)«riS(I) 
10^  FORMAT! 2F1C.^) 

0IN( I )-01h ( I) /lOOO .0 
50  CONTINUE 

C  READ   SOLIC    CtSCPIPfWE    DATA,   M  St)L- IN  1  T  I  Al  POSITILN, 

C  TL-LENCTh,TV— »«10TH,TT-THICKNtSS,Sw-SATu»tTEU  «T., 

C  SEO-SPECifIC    ENEPGr   -lEOUIl-fL    If    INlTW.i:    ELJ«  PAiT 

C  SOL  ID  ,XK-ft€'h    Tfl    SPECIFIC    E^^^CY   COtrf.    AI  SP 

C  VALUfS   ABOVE  SEO 

C  TAU-SOLID    SUfvfACf    C^ffF  Of    SH(  A  ».  ,  F  Hi,  y  -  SL  I  J  1  NG  SOlID 

C  10   PlPr    N/Ll    FRICTins    COi  F  F      ^TAT-S  TM  K    i  iLiu  lu 

C  Pi  PI    HALL    FFiCTICS   C      F  F  . 
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C         soLTiM-Tirr  vnsiTiaM  of  sclio  ih  inflow  profile. note 

C  ZECn  VALUfc    i:  'JIC/T£S    STATIL.^AkY    IN   PIPE    AT  OESICNATEO 

C  SECTION  KSr ! , I NSCL-lsniCATCR  CF  SOLID  VEL OC I TT tCT . 0 

C  INDICATES   THAT  FLT**  ^'ELOCITY  AT   SECTION  NSOL  IS  ASSUrtED 

C  AS  CALCULATION  START  PQIMl, 

REA0I4,102I  HSCL 

IF   INSOL.EO.OIGCTO  7l<i 

REAOMtlOO}  IL«TM,rT.SH 

READM.IO^)  SEO«XK 

REAOf  4«109)Xf AC«XFAC? .T AU«FnOV «FSTAT 

109  F0RrtAT<5F10-^l 

READ!  <i  ,  105  ISOLT  IM,  I  KS  OL 
105         FORKATiFlO.^, 13) 
729  CONTINUE 

KSOLl-NSOL 

IFJINSOL.GT.O  I  NSOL->> 

c 
c 

c  , 
c 
c 
c 

C  CALCULATIOh  OF    INITIAL   CONCITIThS   ALONG   THE  PIPE 

BASED    ON  GPADU^LLY    (VARIED   FLOW   EUUATIONS.  THE 
C  TERMINAL   CTNLITIONS    IN   TH£    PIPE   A^F   dAStu  ON 

C  CRITICAL   DEPTH  Al    «»If»E  EXM. 

CALL   DEPTH«71ME ) 

CALL    PROf  1L<C1H»1I  .H-itHCXN) 

CALL  ASSICNIN) 

IFITlMf.LE.SCiLTIM.ANJ.lNSCL.CT.C)  VSJL>^V(NS0l) 
C  CALCULATION  OP    Tllf   STfP    A^C   LFNoTH  SECTlON. 

CALL  TIMINC(N,VP,CPtVN,CN) 
OX  =  PL/FLO>tT(N) 
DT-DX/(Tf AC«( CN^VH) ) 
DTO=DT 

c 

C  CALCULATION  CF    SOLID    INITIAL  PrSiTIO.-t 

IFINSOLl.CT.C)  XS3L»L)X*FLLAT(hSCLl-l) 
IF<NS0L1.CT.C  )    XPJS  =  lIX«FLCAT(NSCL1-1  ) 

c 
c 
c 
c 

C  OUTPUT   TEST   DESCRIPTION  PLUS    INITIAL  CONIjITIjNS. 

WRITE (  3, 202  ID. RM» SO, PL 

202  F0RMAT(  lHl,/10X,nEST    PIPE    CONF  ICUCAI  I  Qk: , 
1   /lOX, 'OI  ANf  TE  R   «   •»FlO.'i,«    f  .  *  , 

I   iOX, 'MANNING   COEFFICIENT:/-*,?  1C.1t 

1   /10X,«PIFE    SLCFE   -    '.FIO.^,'  PIPE   LENGTH  •  », 

I  FIO.^  ,»  «.•,//) 

IFJHC.GT.HN)   WRITE  n.200)Hf,hC 

IFfHC.LE.bN)    HR  ITE  (  3,  201)»-N,HC 

200  fOR«AT(lOX,*FLOH  SUPERCRITICALf    •.•NO^inAL    DEPTH  -  *«flO.<i» 
I    •    H.«,«    AND   CRITICAL   QEPIh   «    •,F10.'.,'  «.'♦//) 

201  FORMAT!  lOX, 'FLCH   S  Uf<C  R  I  T  I  C  AL  ,    '.•NORinAL   OEPTri  - 

1   F10.^,«    n.    AND   CRITICAL   DEPTH   «    •.FlO.'i,*  «.*,//! 
QiNd  »-QIN<l)*10C0.0 
WR  I  Tf  «  3,  ^03  )tj  IH  <  1  >  tMH  ,CN 
01Ntli-OINU)/1000.0 

203  FORMATdOX, 'INITIAL  FLOW  RATI    -    '.FIO.S,'  L/>.', 

1  •    INITIAL    DEPTH  -    '.FlO.'i,'  N.', 

2  •    INITIAL    hAVt    SPEED   -    '  ,  f  1  0 .    ♦  *  M/S.',/l 
Nl-N*  I 
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WRITE  <3»20«i»L>T,rx,Nl«  TFAC 
20*         fORriAKlOy.'CALCULAriCN   TIME   STEP   -  •.FIO.*,*  S.*. 

1  •LEMCTM   iNCREhLNT   -   •tElO.4,*   «.*»•   HUMBtR   W    HOOES  - 

2  •     TFAC  -•,E5.?,/» 

If  (NSOLI.EO.OI  GOTO  730 
NRITE(3«220) 

220  f ORMATClOXt'SCI.fD  OFSCRIfTIO*  0*TA.«./l 
MftITEI3*221)TL«TM«  T T , S M, NS CL 1 • X SOL 

221  FORMAK  lOX,  •LENGTH  -•,F7.2«*  WIDTH  -'.FT.i, 
7   •    HM.«,«    THICKMESS   -•,F7.2,»  SAlURATtS  HT. 
2  F7.2,«  CM.',//, 

j.  lOX, "SOLID  POSITION,   SECTICN  HUhBER  »»,I3, 
4  •   DISTANCE    FROH  ENTRf    -•,F7.3,*  «.•,/! 

MRITEJ3,222)FSTAT,FMQV,TAIJ 

222  FORHAT{/ICX,«  STATI C  FRICIICH  FACTOR  «  •,F6.3, 

1  •       SLIDING  FRICTION   FACTOR  -  •,F6.3, 

2  /,10X, 'SOLID   SURFACF    SHEA^   COEFF.   '  •,F7.A,/i 
«R1TEC3,2  71)XFAC,XFAC2 

271  F0RHAT(10X,«eCLTAHCY   FACTCR   XFAC  •,F5.2,/ 

I      ,10X, 'DOWNThkLST   FACTOR   XFAC2   '    '  ,  F    .  2  ,  /  ) 
IF(SOLTIM.CT.C.O)   WRI TE ( 3 , 22 3  )  SCLT I M 

223  FORMAT (/lOX, • SOL  10   PJSITION    IN   INFLOn    PROFILE    ^  *, 
I   F6.2,*    SECONDS   FROM  FLOW    I NCC PT I CN . • , / / ) 

TL=TL/1C0C.C 
TH-Tw/1000.0 
TT«TT/1000.C 
SH-SH/1 OOC.C 
730  CONTINUE 

C  OUTPUT   TEST   SirULATIJN  RtStLTS. 

IF (NSOL .GT.O)    GOTO  7i3 
WRITE  (  3,205  )  <  XM  I»  ,1  «  1  ,N1  ,H2) 

205  FO(<NAT(/iey,  '    F0SITT3N   FRf^    •,/,13X,»    ENTRY  N  .  '  ,  fc  x  ,  1  If  7 .  j  ,  /  t 
WRITE  I3,2C*.  )T  INE  ,IHPI  1  »  ,  I-l,Nl,Si» 

206  F0RMAT(2x,  M1^E  -   '.Fb.Z,*    5.',*   D  f  p  r  H  • ,  7  *  ,  •  n=  *  ,  i  i  F  7 . ) 
HRITE{3,2C7)«VP(H,I='l,Nl,N2) 

WRITE ( 3, 3C7)( CP ( 1) , 1= 1 ,N1 ) 
307         FORNAT(  ItJX,  •   FLOW  RATE   L/ S=  '  ,  11  F  7.  ) 
HRITt(3,20E)<CFin,I=l,Ni,N2) 

207  FORMATdtiX,'    VELOCITY      f  /  S -=  •  ,  IIF  7.  A  ) 
203         F0RKATC18X,'    WAVE   SPD.    M/ S  •= '  ,  1  IF  7.  A  ,  /  ) 
C 

uOTO  73A 
733  CONTINUE 

N3-NS0L-2*N2 
KA»=NS0L-1.2 
H5«NS0L*N2 
H6-NS0L-»2«N2 
IF(N3.LE.O)  N3-NS0L-2 
IF(N4.LE.O)  N4,-NS0L-1 
WRITE  I  3,3C^>)  VSOL 
305         FORMATI/^>■«.X,•  SOLia   PG  S  I  T I  CN*  , /5^X,«  SOL  1 0  VELOCITY',/, 
I   5*X,»    -   •  ,F6. 3,»   N/S. • » 

WRITE(3,2  05)XN(l),KN<N3»,XN<N^I,XShL,XS0L,XN(N5»,XNJNbI,*-ifNl) 
WRITE(3,2  06)TIHE,HP(L),HP(h3),hP(N<.),HPUS,HPuS,hPtNS),HP<Nb), 
1  HP(N1 I 

WRlTE(3,207)VP(lJ,VPlN3»,VP(NAJ,VPUS,VPLJL»VPjN5),<'P(Mo>,Vi'(Ni) 
4RlTE(3,30  7lCP(l),aP(S3),GP<N',J,CP5QL,OPSuL,OH(N5),t,P(Nf)).JP(.Nll 

WRITE<3,2C6}CPll),CPtNJ),Cr(N^),C.>'US.CPCi,CP(N'-),CP{'io),cP(Nx> 
73^  CONTINUE 
Z 

c 

C  UPDATE    TINE   ;HC  CHfCK   CALCtLAllLh  LtriiiTn. 

600  CONTINUE 
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501  COHTIHUE 

CALL   TlrlHC(N,VP,Cf ,VN,CNl 
OT-DX/f  TFAC*«CH-»VN1  » 
OTO-DT 

TIME-T IHE*DT 


IFCTIME.CT.TflAX)  COT']  500 

c 
c 
c 
c 
c 

C  SET  Of  CALCttATIOHS  FOR  THIS  Tl^£  STEP. 

OTRAT-DT/DTO 
¥SOL-VPSUL 

IF(DTRAT.LT.0.C6)   GOTO  500 

CALL    iHFLOHfT IHE,QI 

Jl-l 

CALL  INTER^^) 
Jl  =  2 

CALL  LOSS  <VR,VS,HR,HS,M,SR,SS,f ^) 
Jl»3 

CALL  ENTkYITihE) 
CALL  NODAL  I N) 

IFdlKE.GE.SLLTlr.lNj.NSOl  .CT.C)   CALL  SOLlOlTIt'tJ 
IF(TI MF .GT  .SLLT  IM. »NJ . INSLL.CT .( )    CALL    I h S £ R T ( X N , T  I  «L ) 
IF (TI .Of  .SCLT IH. . NSGL .CT . f  I  XPGS=aSul 

IF  (HN.CT.HC  )CALL    E  <  I  T  (  T  I  rtC  ) 
IF{DT .LT.CTC)   GOTO  600 
UT=DTO 

CALL   ASSIChtN  » 

IF(TI       .Gt  .SGLTIM.4'<)  .NSOL.CT  .0)   CfiL-    F  OK  CE  <  F  I  HE  ,  XN  ,F  S  Ol  I 

IF (NSOL.GT  .C)    GCTJ  7il 

Jl»7 
787  CONTINUE 

WRITE(3,206)TlhE,(HP(I),I=l,NI,K^) 

Hkl TE < 3,2  07) ( VP ( 1)  , I  =  I  ) 

HRITE{3,3G7)(CF{I),I^1»N1,>>.2) 

WRITf(3,2fFHCF(I),I"l,Nl,N^) 

GOTO  501 
731  COHTINLIE 

N3»NS0L-2*N? 

N^«NS0L-H? 
,M5-NSGL*n2 

N6-NS0L*2*h2 

IF(N6.GE.hI»  H5-NS0L*! 

IF(K6.CE.NU  ht'=NSaL»2 

IFCNSCL.GE  .H-U  NSOL»0 

IFCN3.LE.C>  h^-NSO'.-l 

IF(N3.LE.0I  h3-NSaL-2 

IF  (NSOL  .E  0  .0)    KRI  T'^(  3t22^) 
22^  FORHAT  I ///30X  , 'SLL  10    HlTHlN  ^    SFCTIOMS    (If    PiPt    D I  SCHAi' G  f  .  *  ,  /  . 

1  30X,'SnLlG    *SStM£3    TJ    CLEAf    P  I  F  I  ,  C  a  L  C  UL  A  I  i  C  N   T  f  «  h  1  N  A  H    ,  S  JL  J I  1 'JN  ' 

2  »/30X  ,  •  C.JH7  IhUES  Win    SPLIT  FRFF    PIPc    Fluh. ',///) 
IFIHS0L,tO.C)WFlI£(3,205)XhJl),)M<Nj),xNiM'i),XPGS»tPjS»Xh(N5)» 

IF  (NSnL.fc  CO    GOIJ  737 
MRITf  J  3,3tc  )  SOLVEL 
303  FOS  HAT  ( >  ,*  SLL  10   ve  L  OC  IT  Y  '  ,/ 5^>,  »    =    •,fb.3»»  H/S,'l 

IF ( VSGL.^T  .C. 0)   GOTO    5  77 


68 


GOTO  ^78 
577  CONTINUE 

WRITE!  3«2C5)XNI  11,XH(  (131  «K»4(f<^  },XPOS»XPOS«  XNlH51,XH(Nbl  «XN(H1  ) 

$78  CONTINUE 

MRITE  f  3«206)T1>1E»HPI  1»  «HP(»43)..f4PIN«  I.HfUS«HPOS«HP  f  NiltHPtNb  ). 
1   HP  INI  I 

WRITEf 3,207)VPt 1I«VP(N3} «VP(N4  >«VPUS«  VPDS . VP! H5 I » V P (H6 ) , VP C N 1 1 

MRlTEi  3t30  7)CP  <ll,ClP(p<3)»QPIHAI  ,CPSOL.OPSOL«QP(hb*  »0PiH6l  t  QP<N1  I 
ttRITE(3«2C8)CPCl)*CP(N3)«CP(N^l«CPUStCPDS«CP(N&)«CP(K&) tCPCNII 
GOTO  501 
300  CONTINUE 
£N0 


SUBROUTINE    TI  H1NC<  N,VP  ,CP«  Vh,CM 

COMHJN/Cr.e/ IhSCl  NSQLI.NSOL 

COHMQ N /CHS/ SOLVE L »HDS .VOS .CDStHPOS , VP  OS, CP  OS 

COMMON /C«1C/X PCS, X  SOL, VP SCL,VSCL ,0P5 JL , HU S , VU S , CU S , HPUS , V PU S , C P uS 
C  THIS    SU3kCUTlNE    IDENTIFIES   THE   HIGHEST    HAVE  SPEED 

C  AND  THE   HIGHEST   AVERAGE    FLCW  VFLlCITY  IN 

C  THE  SlMUL'TtL    FLOW   IN   ORDER   TO  EnSUK;    THAT    THt    T I  ?i  E  STEP 

C  CHOOSFN    IS   THE    SMALLEST,    HfhCE    ENSURING  STABILITY. 

DIMENSION   VF(61),  CPtf)l» 

CN«=0.  0 

DO   1  I»1,N*1 
IFd.EQ.NSCDCCTO  1 
IFICPd  ).CE.CN)  CN»CP(  I) 

1  CONTINUE 
VN-0. 0 

DO   2   I -I.N ♦! 

IF( I.EQ.NSCL)    GOTQ  2 

IF  (  VP  (  I  )  .OF  .VN  I    VN«VP  I  I  I 

2  CONTINUE 
IFICPUS.GT.CN  )  CN^-CPJS 
IF  JCPDS.GT.CM  CN-CP3S 
IF  (VPUS.GT.VN  »  Vh  VJS 
IF(VPDS.GT.VhJ  VN-VPJS 
RETURN 

END 

C 
C 

c 
c 

C 
C 

c 

C 

SUBROUTINE  LOSS(VR,VS,HR,HS,N,SR,SS,<n) 

01  MENS  ION  VF(61),VS<(>i),HR(tll,HS(6It,SR(bl),SS<bll 

COMMON/CMi'/hCHN 

C  THIS   SUBkCUTlNE    CALCULATES    THE    fOUTViLEM    STEADY  STATE 

C  LOSS   BASEL   CN   THE   MA'SNInG  CCIFFICIEnI. 

C 

DO   1  1-2. 

CALL   SHAPHT1ME,HR<II  ,A,T,P£f<» 

SR(  I  )  -  (  (  VM  I  I  «AeSl  VR(  IM«rr*»2»)/{A/P£P.  I«^»1.333 
I  CONTINUE 
HI  «N 
NY-l 


69 


IF(KH.Le.HC) 
IFfHN.LE.HCI  NZ-N*I 
00  2  l-Nt«HZ 

CALL   SKAPE<TIK£,HSf IttAtT.PERI 

SSI  1  )-UVSI  IMABSirS  f  I  >  l*RM**2l>/(A/PERI*»l«333 

2  CONTIHUE 
RETURN 
END 

c 
c 
c 

G 
C 

c 

SUBROUTINE  DEPTHITMEI 
C  THIS  SUBROUTINE   USES   A  SECTION  Cf   ENTkY  TO  CALCULATE  NORMAL 

C  AND  CRITIC/L   FLOW  DEPTHS. 

OIMEMSION   CIh( 30), TIN (30J 

C0MM0N/CM3/NPTS,OIH,r IN 

COMMON/Cm!;/HC  »HN 

CALL  ENTKY(TlflE» 

RETURN 

END 

c 
c 
c 
c 
c 
c 

SUBROUTlNf    IhFLUH(  rue  ,QAV  ) 
C  THIS   SUBRCUTIhE    CALC'JLAltS    iNFLCt-   'ATfS   AT    PIPF    tNTkT  BaScO 

C  ON   THE    ENTFY   FLCH   °^jFlLt    LAT*.    ^.^Tt    IHAI    T  He    0  LALCoLATtJ 

C  IS   AN   AVEfrAGt    VALUE   rOR    THIS    Tift  STtP. 

riMENSION   0(  30  )  ,T  (  1^^  I  ,QX  t  jr  ) 

COMMJN/C^l/r  T0,L1  ,0<,D,SO,PL,SFO,)"r,-<<^ 

C0MM0N/CM3/NPTS»O, T 

Tl-TIMF 

TO«=TI  ME-DT 

J=I 

TX  =  T1 

DO  3  I-1,NFTS-1 

IF  CTX.Ct  .T  (  I  )  .  ASC.  K.LT  .T  t  Ml  I  J  LQT:^ 

3  CONTINUE 

4  0X(  Jl  «C(  I  )  ■»<&(  1-HJ  -0(  I»  )»(  TX-T  <I)l/(T(I«l!-ni)) 
QAV«ax(lJ 

P.ETURN 
END 

C 

c 
c 
c 
c 
c 
c 

SUBROUTINE    SHAPE (T I1£  ,H,A,T,PEP > 

C 
C 

DIMENSION    OIN  no  )  ,  T  lr<  (  30  > 
CuMHON/Cfil/C'?0,LT.l)X  ,U,S0,PL»5I  (.  .XK.K.i 
COMM0N/Cfl3/NFTS,O[N,TIN 
C0M»<JM/Crt7/CL 

C  THIS    SUp.P.CUTIht    CALCJLATIS    fLfW   AC"  ^  t ,  S  U^;  f  A  C  £    wILiTH  ANLi 

C  HETFED   PL^If-LTtR  BAScU  FLlW  tcPTH 
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C  THIS  SUBROUTINE   ALSO  C»LCU1.A1E5   MATER  SURFACE  fPOflLE 

C  AMD  SETS  UP    PASE   CJ><aiTIOMS    ALONG  THE   PIPE    AT   I  I  Hfc  ZERO. 

R-0/2,0 

|F«M.LT.R)T»^f  TA-2.  0»A T  ANtSCRI  IH«<0-H*  I  /IR-HII 
IFCH.EQ.RITHETA-Pl 

iF<H.CT.r.  IT^E  TA-PI  «2.  0  «  ATAM  (  (H-Ri/ (  S3  R  T  C  H*  (D-H))!  ) 

A-( (04*2 ) /B.O)*(lHETA-SINCTh£TA)  J 

fER-D*THETA/2.0 

1-2. 0*<  IH«<{>-HI  »♦♦<). i  I 

1F«TIME.GT,0.0»  COTO  1 

0-OIN(  1  > 

C0H-RH»«2/S0 

HCRIT-1.0-<0»»2l»T/f i»A**3  I 

HNnRh-1.0-|0»«2J»C'lN/  ((A««3.333»/IPEil*»1.33  3n 
DL-HC  R  1 1/  t  UNO  RK.*SO  } 
1  CONTINUE 
RETURN 
END 

C 
C 
C 
C 
C 
C 

SUBROVTINt  »<AVSPD{H,C) 

C 

C  THIS    SUBKCUTIhE  CALCULAT£5    WAVE    SPFEO    BASED  UN  DEPTH  AN'j  CROSS 

C  SECTION  ShAPE. 

CALL    SHAPE «TlHE,h, AREA, T,PERJ 

C=SQRT(9. fcl*AREA/T) 

RETURN 

£ND 

C 
C 

c 
c 

SUBROUTINE  PROFILCO,HN»HC,XN) 
C  THIS    SUBKCUTIHE    CALCJLATES    ThF    INITIAL   WATER  SURFACE 

C  PROFILE    BASED   ON  CRITICAL    EEPTH  AT   PIPE  cXiT. 

DIMENSION  HPteD.VPCiD^QPCtlNCPCfal) 

DIMENSION    X (30 » tOF" ( JO » ,XN (61 ) 

DIHENSTOn   X 1( 3C  I  tOE* I ( 30J 

COMMON/Cnl/Cl  0,DT  ,3X,0,S0,PL,SEC,XK,k:< 

COMMCN/Cr.i  /'N  ,  TPIAX 

C0WM0N/CMt/OP,VP,HP,CP 

C0MMQN/Ch7/tL 

COHnON/CMJi/  INSOL  ,H  SOi.  1  ,NSCL 

C0MMUN/Cri9/SCLVEL,HDS»VDS,CrS,^'PLS,VPi)S,l-P0S 

COMM0N/C«10/XP0S,XSQL  »  VPSt  L»VSOLtOPSjL,hUS,VUS,CUS,HPL.S,Vr'US.CPUS 

c 

IFCHN.LE.HCJ   COTO  900 

DH-(HN-HCI/30.0 

IS-1 

Hl-HC 

CALL   SHAPE (TIME, HI, A, T, PER ) 
X(  ll-PL 
OEPd  »-HC 
SL-0. 0 

C    '         .<ATER    SURFACE    PROFILE    C  AlCULAT  I  CnS  . 
00   80  I-l,;'CC,2 
1S-IS»1 
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M2-HC»DH«FL0A7  U*ll 

H3-KC*0H»FL0AT(XI 

CALL  SHAFE(Tint»Hl«A,T«PEtl 

OLl-OL 

CALL   SHAPE  (TlflE»H2«A,T*r€«t> 

DL2-0L 

CALL  SHAPEfTirE*H3«A» TtPERt 
0L3-0L 

OXP-OH*tOLl«DL2«^.0«JL3 1/3.0 

SL-SL-DXP 

«l-H2 

IFISL.GE.PLI   CCTO  91 

XUSI-PL-Sl 

IF(Hl.CE.HH)   GOTO  83 

DEPflSI-Hl 
60  CONTINUE 
81  XIISi'0.0 

NIS-I S 

IF  CHI  .  CE.»iN)   GOTO  93 
OEP( IS l-Hl 
GOTO  8^ 
83  OEPtIS)-HN 
GOTO 

8A  IF (X( I S).CT.C.C1   C3TJ  65 

GOTO  6fc 
85  X(1S)=0.0 

DEP(iS)»riN 
■  NIS«I S 

66  CONTINUE 

e 
c 

C  INTERPOLAT  ICh   PEOUIREO   TO  CFTtPMNf  0£PTH  AT   EACn   NQjt  . 

900  CONTINUE 
DX»^PL/FLOA  T<h  ) 
XN(l)-0.0 

00   87  I=2,N-»1 

XN{ I )-XN(I-ll*DX 
87  CONTINUE 

IF(HH.LE.HC)   GOTC  901 
S2=NIS»1 

DO   9*1  J-1,N1S 

N2-N2-1 

Xl( J»=X(N?J 

OEPK  Jl-0EP»N2  ) 

MRI  T£  I  3t93  >X1  (  J  ),0E  PI  (  J  » 
93  FORMA T{ lOX  ,2f 10.51 

9^  CONTINUE 

03  *»5  J"1,NIS 
XOI-XKJI 
OEPC J)=0£P1 ( J» 
95  CONTINUE 

HPCD-OEPtn 

HP(N*1 l-HC 

00   98    1-2, N 

DO    89  K-1,HIS-1 

IF  (  XN  (  n  .C  T.X  tK  )  .Ar«n.  Xr<(  1  )  .Lt  .y  <K*  U  )    CGI  (J  9o 

89  CONTINUE 

90  HP(n-OEP<»<)«IDEP(K»l)-DtP(K))#(X><fI)-X(^l  »/lX<K«l)-x(K)) 
83  CONTINUE 

C  SET   UP    PASf    CLhCITION*    At    Tll-£    2 1  iiD . 

901  CONTINUE 
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CALL  :>MAP  t  I  ;  1  r  t  ,Mr  1 1  5  »  A  »  I  »  FtR  1 
¥Pt I  I -C/A 

QP  t  n  -Ci«ic 00. 0 

CALL  MAVSF0(1!P(I)(CPI  !t) 

CONTINUE 

IF  INSOLl.EQ.O  »  COTO  99 
HPUS- HPINSCLI » 
MPDS='HP<NSCL1J 
VPOS" vr<N5CLl » 
VPDS-VPIH5C1L1) 
CP'JS-CrtNSLLl) 
CPOS-CPCNSOLII 
i2P50L-OP<HSOLl) 
99  CONTINUE 
RETURN 
END 

C 
C 

c 
c 
c 
c 
c 
c 

SUBROUTINE  INTER(N» 

C 

01  MENS  ION  V<tl),H<61),C<f>l),VK(tl),H<(6l),CR(61),VS(on 
DIMENSION   HS(6l),CS(tjl),SF  <til,5S(fcH,XN{bl) 
DIMENSION    >Sl.6l)  »XPf  bl) 

COMMON/Cni/CTO,DT,Dt,0,SO,PL,SFC,XK,<n 
COHM0N/Cfl^/V,H,C,VP,HR,CK,)(l-,SP,VStHi,C5,XS,SS,XN 
C0MM0N/Cn5/HC  -thti 
CONMON/CMb/  INSOL ,N  SOL  1 »NSCL 

COMMON/CM^/ SO L VEL » MOS , VOS , CCS ,HPLS,  VP  OS ,C  PDS 

COMMON/CNlO/XFCS,XSnL»VPSCL,VS(L,OPSaL,HL,S,VuS,CuS,HPuS,V->US,CP'Ji 
C  THIS   SUBRCUTINE   SETS   UP,    PY   i  NT  t  k  P  OLA  T  I  L  N  ,    Trtt    BASE  ^:ONOiriJNi 

C  FOR    THE   NEXT   TlrE  STEP. 

C 

THETA=DT/DX 
N1»=N*1 

IF(  VSOL.GT.O.OI   CCJ  .1  6 
COTO  7 

6  XA*=XSOL-OX 
XC-XSOL 

DXA«XA-XH( NSOL-1 ) 

VA-V( HSOL-1 » ♦ ( VUS-V (S  SOL-1 1  l*D^  */OV 
HA»H(NSOL-U-»  (HUS-H(NS0L-1  )  )•  CXA/Or 
CA-:C(  NSOL-1  !♦  (CL'S-C(NS0L-1  )  >*L>*/OX 

VCU«  V  (NSOL-l)*(VLlS-VJ  NSOL-  1  )  )  «  (  >  Si'lL-XH  «  NS  OL-i  J  )*/  U  K 
MCU-H  I  NSDL-1  )  ■»  (HUS-H(  NSOL-in*(XSnL-AN<NSuL-l))/LX 
CCU-C  (  NSOL-D*  (CUS-CINS0L-1))«  (  >  SOL-XN  (  NS  JL-i  i  i  /  u  K 
VCO-VDS*( V (HS0L*1 » -V3S  >  *( XSLL-XK(NSLV  ) ) /tX 
HCO-HDS*(HtNSOL-»l»-H0S)»(X50L->N(NS0L  »  )/Ua 
CCD-CDS*  (C  t  NSCL^  ll-COS  »  •  (  X  SCL -X  ^  (N  S  OL  M/LX 

xe-xsoL-»ax 

0XB-XB-XN<NS0L-»1I 

VB-V(NSOL«l>*iV(NSOL*2)-VC^SUL<])»»D<:i3/D* 
KB'=H<NSnL*l)«(h(NS1L»2)-H<NSL;L*lJ>^OX6/DX 
CB-C<NS0L-»H-»<C<NSlL*2l-C<hS0L-»i)>*DXo/0X 

7  CONTINUE 

00   1  i-?«Nl 

IFt  I.E  O.NSCL./.M0.VSO;_  .CT.C  .C)    Ctl"  6 
tOTO  9 
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8  ¥<II-VCU 
H(  n-HCU 
Ct  I  l-CCU 
»C I-ll-VA 
MCI-l >-HA 
CU-1  l-CA 

9  COKTIHUE 
IF(VSOL.LE.O,0.AH0.I.EO.NSOLI    V(  IMVUS 
IF«VSOL.LE.O.O.AND.I.EQ.NSLL>   C«  n=CUS 
IFIVSOL.LE.C.0.AN0.I.EC.N5CL)   h( I)=HJS 

If  «V50L.L£.0.C.#N0.I.EQ.r*S0L*l)  VCI-1»-VDS 
IFJVSOL.LE.O.O.ANO.l .EQ.MSOL*! »  C(T-l»*CDS 
IFI  VSOL.LE  .C.O.ANO-  I  .  EQ.MSOL*  1 J  fit  I -1 1- HO  S 

vR(n-=cvii)*ThE7A<'jc{n«v(i-i)-vui*cti-i)n 

I  /C  l.O^TME  TA*1  VU  I-i  )  ♦€  <  n-C<  l-l)  )  ) 

CR(  I  )  -  (Ci  1  )*<  1  .C-VR(  I  )«THET  AMC  (  1-1  I*  VRt  I  1  »IH£  TA) 

1   /(  1,0-»C(  1  )*T^  ETA-C  (  I-l  )  »rHE  TAJ 

HRC  I  l-HU  )-<h(l  l-H(  I-l  H»1^ETA«<  VS  (1  »»CR<  1  )1 

1  CONTINUE 
IF(HN.LE.HC>    CGTO  ^ 

00  2  I^ltN 

IFd.EO.NSCL.AhD.VSOL.GT.O.OJ    GOTO  10 
GOTO  11 

10  VdlxVCD 
H(  I  l-HCD 
C<I)-CCD 

vn*ii-v6 

M(I*lJ-Hs 

c  f  1*1  )'=ce 

11  CONTINUE 

IFCVSOL.LE  .C. C.ANO. I . E3,NSCL)  Vtl  }=V3S 
IF  I  VSOL.LE  .CO.  AND.  I.EQ.NSCLJ  Cn)=Cjj 

1  FC  VSOL  .LE  .C.  0.  AND.  I  .  EQ.NSCL)    ^'(  H  =HjS 
IF(VSOL.Lt.C.O.AN;m.£Q.NSCL-l)    V(  I<-..»=VuS 
IF(  VS  OL  .LE  .  0.  O.ANO.  I  .  EQ.N5,LL-1  )   C  (  I  +1  ►=Cb  5 
IFtVSOL.LE.CC.AhT.I.  E  S.-sSCL-l  )    h<  I      i  =Ht  S 
VS  m  *  (  V(  I  )-Th£TA*  (  V(  i  )«C  (  1*1  )-C  (  H»»  «  1*1  )  »  I 

I   / r 1.0-THET A«( V  I  I  )- V( I  ♦! t-c ( n (1 ♦  1) J  » 

cs<n»{Ccn-»vs(n>THETa«({.<i)-tu*in> 

I   /(1.0*7HETA«(C(n-C(l*l)  »  ) 

HSUJ=H(U-»ThETA«(*<S(i)-CS(I)J*th([)-H(I*l)) 

2  CONTINUE 
GOTO  3 

*  DO  5  l'2^h*l 

IF  (  VSOL  .CT  .  C.O.ANO.  I  .  E  Q.NSCL*  1  »    CvlT.,  b 

IF(NSnL.GT.C,0.AN0.!.EQ.N5CL-»l)   CTU  20 

GOTO  90 
20  V(I-1}«VDS 

M« I-l » -HQS 

CC I-l  )-COS 
80  CONTINUE 

VS(I)'lvnj»l  1  .0»TH^TA«C(  I-IJ  J-V(  I-1)»THlTA»C(1)I/ 

I  <  i.o*thet*  ♦  J  V  (  n-v  ( i-ij*c  u-1  »-c  ( 1  n ) 

CS(I)  =  (Ctl)»VS<I)»(C(IJ-C<I-l)))/<l.J»TH,  TA«(CU-il-t<I))» 
HS(H-H(I)-<hCn-H(I-l)J»THfT^«»yS(I  )-CS  (  I  )  ) 
5  CONTINUE 

3  CONTINUE 

If ( NSCL.EC.OI    CCTQ  13 
IF (HN.CE . HC 8    GOTO  18 
DO    17  I-l,hSLL-l 
IF  (  VU  I  .01  .CU  )  )   CHTl  17 

vs<i)-(vin-in£TA»(v(i»»cn-»i)-c<M»vii*i))) 

1    /(1.0-UvT**(V(I)-V(i  ♦ll-C<l)*f  <iM  >>l 


74 


CSf  II-CC<n-»VS<l  I*TM^TA»ICI1|-C«l»ll>l 
/(1.0«THC  TA«(C(  I  l-r.  I]  ♦!  )  i  I 

HS<I  l-^Hf  I  l*ThETA*(  VSt  I  >-CS(I>  M(H(  II-HCI«il> 

CONTINUE 
CONTI NUE 
RETURN 
ENO 


SUBROUTINE  ENTRYtTIHEl 

THIS   SUBROUTINE    CALCJLATES    THE   UPSTREAM  BOUNDARY  CONDITIONS 
AT   EACH  TI^E    STEP   BASED   ON  A   KrChN    MfLOW  PROFILE. 
01  HENS  I  ON  0P( 61 ) tOINi  30) «TIN4  3C) tVP(b I) ,HP (61 ) tCP I 61i 
OINENSION  V(bH,H(6n 

DIHENSION  VR(Cl)«HR(olltCP(&ll«VS(61itHS(6l)tCS(6l)tSS(bl) 
OINENSION  5PC61),X«(*>n,XSCtl),XS<'>ll,C(t.l) 
C0MM0N/CM1/C7  0,DT,0<,  D,SO,PL,  SEC  .XK,<n 
COMMON/Cfli  /Nt  T^AX 
CONHON/C«3/NPT5»CIN,r IN 

COMHON/CM^/V,h,C  »VS»^'<  ,CR,yf<  »SP  »  VS  ,HS  ,CS,XS,Sb»XN 
C0NMQN/CM5./hC  ,HN 
C0MM0N/C«t/Ot ,VP,^P,CP 

C-9.81 

C0N*RN»«2/SC 

|F(T1MF.CT.C.0)CALL    INFLOWdl^F  ,C) 
IFCTIME.CT.C.C)    CQTl  600 
IKTIME.EC.O.0)O-3fN(  I) 
CALCULATION   Lf    CRITICAL  DEPTH. 
UP«D 
0N=0.  0 

HC*(UP*DN) /2.0 
CONTI NUE 

CALL   SHAPE »TINE,HC,ARE A,T,PER) 

HCRIT-l.C-(0**2)«r/rG«AR£A««3» 

IFIHCR  ITia^'^.S 

DN-HC 

GOTO  6 

UP"HC 

MCN-(UP*1)N  )/2.0 

IF  (ABS  (  (HCh-HC  )/HC  I  .L  t  .0.  CCD    CCTD  8 

HC-HCN 

GOTO  7 

MC'HCN 

CONTINUE 

CALCULATION  LF   NORMAL  DEPTH. 

UP-D 

ON-0. 0 

HN-»UP*DN 1/2. 0 
CONTINUE 

CALL  SHAPE«TlME,HN,A8EA,T»PEk) 

HNOPM  -  l.O-l C««2>«C0S/ ((AKtA»*3.353>/(Ptk»*1.333>) 

!F(HNORH)  10«lltl2 

Dh-hn 

GOTO  13 

UP  "MM 

HNN«( UP •UN  )  /2  .0 

I  f  (      ^  (  <  H-^*  -t     »  /  »  s  »  .1  f  .  ^  .    r  1  )    m  i    J  <, 
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HN-HNN 
GOTO  9 
HH-HNN 
11  COHTIMUE 

IFf TIME. EO. 0.01  COTQ  900 

C 

c 

c 

C  CALCULATIOW  OF  BOUNDARY  OEfTM. 

600  CONTINUE 

IFtHN.LE.HC)  GOTO  700 

UP-0 

DM-CO 

HB-(UP*DN) /2.0 

15  CONTINUE 

CALL    SHAPE ITIHE.H3, AREA, T»PER > 
X3-G/CS< 1) 

X*-VS{1I-G»(S5(1)-S0) »DT-X3*hSfl) 
HFL0H-O-IX*i«X3«HB>»A»^E  A 
IFIHFLOKJ It, 17,16 

16  UP=H3 
GOTO  19 

18  0N«H8 

l<5  HBB«(  UP^Df*  )/2  .0 

IF { ABS( <HEe-hB  )/H^ » .LE .O.CCl  )   CCTD  ZJ 

rtS'HQ 3 

GOTO  15 
20  HB»HbB 

17  COHTiNUf 
HP(1)»=HB 

VPCl)  =  X',*>3*hP(ll 

CALL    SHAP  t  (TI  hE  ^H"?  .ARE  A  «T,  P£K  ) 

OP'1}=AR£A«VP(1»»100j.O 

CP( 1)  =  50kT (C«AREA/T  > 
500  CONTINUE 

GOTO  800 
700  CONTINUE 

C  CALCULATICN   CF   NCP1A_  DEPTH. 

UP  =  P 
DN=0.0 

HB«=(  UP-»DN  )  /2.  0 
7<5  COiUINUfc 

CALL    SHAPE (TI ^E,h9,ft^EA,T, PEP) 

HN0RM'=1.0-(C««21«C3N/((AkfcA*«3.3  33  >/(PEk«»1.333)) 
IFIHNOPrOfcO,  1,82 
30  ON»HB 

GOTO  83 

82  UP-HB 

83  HBB»(UP*ON  .0 

lF(AaS( (riEf-he)/H3 I.LE.0.001)    OCT :  bt 

HB«H6B 

GOTO  79 
8<i  MB-H6B 
81  CONTINUE 

HP( l)«H8 

CALL    SHAPE  <T  IME  ,l-!3  ,  A.^E  A,I,  P£k  ) 
VP ( 1» -Q/AP  f  A 
Cf  (i  }-Q«iCCi'.0 
CP(1»-S0RT«C»ARE*/T» 

800  CONTINUE 
RE  TURN 
END 

C 
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c 
c 
c 

SUBROUTINE  NCDALIHI 

OIKENSIC<  OP(6l}tQIHI  30ltTIHI30l«VP((>ll»hPf6ll,CPf 611 
OIHENSION  Vf6l}tH(6M 

01  HENS  I  ON  VK(611«HR{61>«CR(61 )tVSf61)tHS(bl)tCS(611«SS(61) 

DiHENSifi  !;r(6U,x?(si  ),xS(6i>«yNCM>tCft>ii 

CCMKGN/CK1/D10«DT,0X«D«  S3»PL tSFC*XK«Rn 
CaNnON/CN3/hPTS«OI Ntf IN 

C0MM0N/CR4/V,H«C«VR,^1R«CR*XR«SP«VS*HS  ,CS*XS  «SS«XM 
C0«M0N/CH5/hC  »HN 
C0HH0N/CN6/OP , VP,HP ,C  P 
C0HM0N/CH8/lNSCL,NSOH.NSOL 

C0M«0N/C«10/XP0S,XSai.,  VPSOL»VS0L,OPSaL  fHUS»VUS,CL»i,HPUS»V  PUS,  CPUS 
C  THIS   SUBROUTINE   CALCJLATES    THE    FLOW  VELOCITY  AND  OEPIH 

C  AND  HAVE    SPEED  AT    EACH   OF    THE    ►CCES   iETMEEN    IHE  uPSTkEAI 

C  AND  OOHNSTREAfl  BOUNOAklES    BY  SCLUriOi   OF    THE   IWO  HAVl 

C  EQUATIONS. 

C-9.81 

N2-N 

IF(HH.LE.HC»  N2-N*l 
00   1  I>2,NZ 
IF<I.EO.NS0L)    GOTO  I 

IF(VSOL.GT.C.0.AND.I.Ea.NSOL-m   tOTO  1 
X1-G/CR< I  ) 
X3-C/CS{i  » 

X2-X1*HR(  I  MVRCI  >-G*J  SR(  I  )-S0  )*CT 
X^  =  VS { 1 l-X3«HS{ 1 )-:*( SS I  I l-SC)«DT 
HP(n  =  (x2-XA)/(Xl»XlJ 
VP(I  )=«X1*X3*hP  (  I  ) 
CALL   HAVSPDihP  11  )i C»(  I) I 
CALL   SHAPf  <TIrE,H5>  (I  »  ,X1,X3,X^) 
OP(I  I-X1«VP<  I  >«1000..} 
I  CONTINUE 
RETURN 
END 


C 
C 
C 

SUBROUTINE  FXITITIHE) 

01  MENS  ION   OP  Ifcl)  ,yo  {(s  I  )  ,HF  (fcl)  ,CP(  Ml 
DIMENSION   VJfcll  ,H'((i>ir  .     •-  • 

OIHENSION    VRl61),H<^<bl>,CR(61  ),VS(Sl)  ,HStoH,CS(cl),SS(6l) 
DIMENSION   SPI 61)  «X^(6I  )  ,XS<fcl 1 ,>N( 61 ) tCIbi ) 
COMMON /C M 1 /DT C , DT, D X, 0, SO, PL, SfO.XK,R« 
COMMON/CM?/N,THAX 
C0MM0N/CM5/HC,hN 

COlMON/Ch«t/V,h,C,VR,HR,CR,>»i,SP,VS,HS,CS,XS,SS,XN 
COMMON/CMfc/OP  ,VP,Hf»  .CP 

C 
C 

C  THIS   SUBROUTINE   CALCULATES   Tht    FLO«   DEPTH   AT   PIPF  OISCHARCf 

C  BASED  ON   THE    ASSUrtfll^N   OF    CRITICAL    Dt'PTu  AT    SuCh   A  BUUNOakY. 

C 
C 

G-9.81 

I 
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J<1-G/CR<H-»1I 

X2-Xl*HR«  K*IMVR(H»1I-C«CSR««*1I-S0I»0T 

UP«0 

ON-0.0 

HB-CUP«0NI/2.0 

CALL  SKAPE(TinE«K3«*,T«PERl 

HEXII-1.0-(AeStX2-Xl«HB)9CX2-Xl*H8l )*T/(C*A) 

JFfHEXITll«2t3 

OH«KB 

COTO  4 

OP  -HB 

GOTO  4 

HB9-I UP*0H)/2 .0 

IF(ABS((HCe-MB)/H3I.LT. 0.00001)  COTO  2 
H6-HBB 
COTO  5 
CONTI HUE 

VP(H* 1 J-X2-X1 ♦HE 
MPfM*ll»He 

CALL   SHAPt<7IHE,H9,A,T,PEF!  ) 

CALL   WAySPr<HB»CP(N  1»J 

OP<N»l»»100C.O«A«VP(t^-»l) 

RETURN 

END 


SUBROUTINE  SCLID(TIHc» 


01  KEN  5  I  ON  CP(fcl)t3lN{30»,TlM3ri,VP(ol»,nP(6l)tCP(t)li 
DIMENSION  V(tl),^{6n 

DI  MENS  ION  VHtll,H?J-,l),CF(tll,VS{ftl»,.HS(&l),CSIfcl),iS(blJ 
OlflENSION  FS(2),S'<(*>l),X,'»fcl),XS('SU,XN(ol),C(tl) 
COri^'ON/Ch  1/C7  0»Dl,D«',D»SO,FL»SfC,>''»'(i 
COHMOh/Cni /h ,T^AX 

comhon/cm3/npts,cm,t  in 

COHMON/Cf"./V,h,C,VR,HR,C(<tXk.St,VS.HStCS,Xi,SS,XN 

COI1HON/Cn5  /HC  thN 
COMMON/CMt/CH  ,VP,SP,-P 
COMMCN/CMt/ INSCL,HSDL l,NStL 

COMMON/CMS/ SC L VELtHOS tVOS *  CDS* KPLS , VP  OS , CP  05 

COHMON/CrilO/>  POS  t<>  Ol  »VPSOL»VSri  tOPSOL.MUS  »VUS,CU  S»HPUS  t  Vi>US»CPUS 

CO>^haN/Cni  l/TL,Tt<t  TT,  5  H 
£  THIS    SUBkCUTIhE    CALC'JLATES    THE    fLOw   DEPTrt,    LEAKAGE    RATE    ft  ND 

;  KAVE    SPEED   AT    THE    lOi^NSTKtAt   BDUKPiPy    FOkMED   BY    AN  INITULLT 

C    '  STATIONARY  StiLlC. 

C-<>.8  I 

C0N«Rh««2 /5C 

QPSOL-OPSCL/1000.0 

IZ-0 

SEV-SEO 

4-1 

HO-HUS 
6C  CONTINUE 
J-  I 

aoOO       CALL    SHAPE  <TIhE,HO,Al,T,P£»;» 

X1-C/CR(NSCL) 

X2-XI *Hk  <  hSOL  )«VP( HSJL l-C«<SR (hSOL  »-S0 )»ul 

SE«HJS*(VL'S'>«?)/t.'.0>G) 

IF«SE.LT.SEO>   GOTO  fc 


o5  tun  I  I MUt 

U-(X2-VP5CL)»A1/XIC 
B-X1»A1 /XK 

r-«X2»*2»/l2.0'»CI-SE>-VPS0L«X2/C«IVrS0L*»ZI/C2.0»CI 
Z-Xl»VPS0L/C«i.O-X2»Xl/6 

Z3-2«0<  W*7 
Z2-Z»»2*2-©*l»»y 

zi«2.o*ir*z«e 

ZC-Y«ABS<  T>-U 

f-Z^#H0*#4*Z3*H0*»3*Z2*HO««2*Zl«HO*Z0 
FS<J»-F 

IFIJ.EQ.21  CCTO  7000 
J-2 

H0-H0*0.0001*M0 
GOTO  BOOO 
7000       t>f-CFSC2  »-FS«  lU/t^O/  1 0000. 01 
OH-FSf 1 l/DF 
Hl-HO-OH 

IFIABS  MH1-HC)/H0)  .LE  .0.005)   CCTO  70 

IFCH1.GT.D)M1-1.1*H0 

HO-Hl 

iz-rz*i 

GOTO  60 

70  CALL    SHAPE (TI NE,H1,42,T,PEF ) 

IF(ABS(<A?-A1)/AU  .Lr. 0.005)   CCTG  80 

Al«(Al<^A2)/2.C 

H0*Hi 

GOTO  65 
80  CONTINUE 

MPUS=H1 

VPUS»X2-Xl«m 

QPS0L»VPUS*/2 

IFtVPUS.LT.VFSOL)   GOTO  37 

GOTO  5 
37  CONTINUE 

VPUS» VPSOL 

HPUS-(X2-VPUS )/Xl 

CALL   SHAPEfTIhEtHPU^, A2TT,P£k) 

OPSOL«VPUS*A2 

GOTO  5 
5  VPUS=»VPSOL 

HPU5-  IX2-VPUS  )/XI 

CALL    SHAPE<TlH£,HPUS,A2»r,PEK ) 
QPSOL- VPUS*A^ 
GOTO  5 

5  CALL   HAVSPD<hPUS,CPUS ) 

SEO-SEV 

OT-DTO 

GOTO  006 
•«06  COHTINUE 

IFIHN.LE.HC)   GOTO  dOal 

GOTO  8062 
8061  CONTINUE 

Q-OPSOL 

IF«O.LF.O.C)    0-0.05/1  10 . 0*  10. 0*»-M 
C  CALCULATiCH   OF   hOKflAL  DEPTH. 

UP-O 
ON-0.  0 

HB-(UP»DN)/2.0 
7«J  CONJIHUL 

CALL   SHAPE  lTlr<E,HB»»K£A,r,Ptk» 
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MN0Rrt-1.0-IC*«2»*C0M/(CARE**«3.333>/tPER»»1.333ll 

IF  (HNnRH}<)0«91«92 

90  08f8B93 

92  UP-Hfl 

93  HB8«>(UP«0H)/2.0 

IFCABSt (HeB-HB}/H3»  U.E. 0.001)   GOTO  <M 

H8"HAB 

GOTO  7«> 

94  M8*M8B 

91  CONTINUE 
HPOS-HB 

CAL*.  SHAPE(TlhE.Ma,A«EA,T,P£RI 
VPOS-O/AREA 
OPDS»0*1000.0 
OPSOL-OPOS 
CPDS» SORT {C*AREA/r I 
GOTO  934 
8062  X3=G/CSCNS0L) 

X4»VS«  NSOL  l-&«ISS<  HS  JL  l-SO  l*DT-)(  3*HS<  NSOL  I 

IFIOPSOL.EC.C.C)   C:T3  932 

UP-D 

ON-0.0 

HB-CUP*ON ) /2. 0 

15  CONTINUE 

CALL    SHAPt  (TIME  ,M3,AREA,T-.PEk  » 
HFLOH*QPSOL-( X4»X3»M8  )«AKFA 
IF(HFLOMI lt»17,ie 

16  UP^HB 
GOTO  19 

18  ON=HB 

19  KBD-^t  UP-»DK  )/2  .C 

IF CAe5( (HBE-HE)/h3l._E.0.C01»    CCIU  ZO 

MB«HBB 

GOTO  15 

20  H8=H53 

17  CONTINUE 
HPOS^HB 

VPDS»=XA*X3*HP  LS 

CALL    SHAPE  !TIf«E,h3, AREA, T^PEkJ 
QPSOL -ARtA«VPDS«10':  .0 
CPDS*  SQPT<  G»Af<EA/T  » 
GOTO  933 

932  CONTINUE 
riP0S»-X4/ X3 
VPOS-0.0 

CALL    HAVSPDOiPCStCPOi  ) 

933  CONTINUE 

IF(VPSOL.CT .C.OI  GOTJ   931  • 
GOTO  935 
93*  CONTINUE 

OXX-DX/< XSLL-XhCHSOL-1 n 

VPINSOL  J»VPINSCL-l » ♦ J  X X« ( V  PUS-VP < NSOL - 1 U 
KPtNSOL }-hP<hSCL-l  ) ♦ J  XX ♦ ( HP US-HP < NSOl -1 » I 
CP<NSOL  » -CP  <NSCL-l »*OXX»iCPUS-CP<N>OL-ll) 
0XX«»XN(NSCL«l)-XS0LI/(XN<r<SUL»2)-XSJL» 
*PJNS0L*l»»VPCS*DX<»(VP(NSCL«2)-VPnS) 
MP  (NS0L*1  J-HPL)S*DXX»(  HP(HSCL«2  l-HPHSi 
CP(N>0L«])»CPDS*&X<»iCP<HSLL'2>-CP0SJ 
935  CONTINUE 
Rf  TURN 
END 

C 

c 
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SUBROUTINE    INStRTI tM, T  iMf I 

DIMENSION   XHfbl  )  ,aP(M)  ,Vf  (61|,Hf  C611«CPCbll 
C0MMQN/Ct16/0P  ,VP,Hf  ,CP 
CONMON/CHfc/ INS0L,NSDL 1 «NSOL 

C0MM0N/C«<^/SDLV£L»ri3S  .VDS.CDS,HPDS,VP OS»CP0S 

C0NK0N/CnlO/XPOS«XSaL*VPSOL*VSCL«OPS0L«HUS«VUS«CUS«HPUS«VPUStCPUS 

c 

INSOL-0 
HSOL-NSOLl 
VPUS- VPIhSCL ) 
HPUS-HP( NSOL) 
CPUS- CPCNSCL ) 
*PDS« VP(NSOL) 
HPDS«HP(NSCL ) 
CPOS«CP<NSOL  i 
OPS0L-0P(lsSCL  > 
VSOL-VPUS 
VPSOL-VSOL 

OX-XNINSGL  »-XN<HSOL-l » 
DXX=DX/4  XSCL-XN(HS3L-1  )  ) 

VP (NSOL I -VP (NSCL-l I ♦axx*( VPUS-VP4NS0L-II ) 
HP(NSOL)'HP  <hSCL-l  >  ♦l)XX*  (hPUS-HP  (NSOl-I)  ) 
CP  (NSOL  l-CP  <NS0L-1  I  ♦JXX»(CPtjS-CP<NSOL-i  I  » 
DXX=(XN(NS0L-H)-XSaL»/(XN(NSLL*2J-)(SuL> 
VP(HS0L*i)«VPLS*Dx:X*lVP(MSrL*2  )-VPOS» 
HP (NS0L*1) "HPUS^CX  <«( HP( NStL*2  )-HPDS) 
CP(NS0L*1  )«CPDS»DX<*(  CP(NSCL-»2)-CP0Sl 
RE  TUR  N 
END 

C 
C 

SUBROUTINE    f 0RCE(TI1r,XN,f SDL) 

C 

C  THIS   SUSRGUTINF    CALCULATES   ThE    fCCE   ACTING  ON   THE  SOLID 

C  AND   DETtPhlHtS    ITS   VELOCITY   AND   POSITION  AT   THE    tNO  Uf    1 n i 

C  NEXT   TIME  STEP. 

C 

OIMFNSIOS  XN(fcl) 

COhMON/Crtl/DTO,DT,0)(,0,SO,PL»SECtXK,Rn 
CrjM-0N/Cn2  /N  ,  ThAX 
COMHON/Cie / INSCLtN SOL  I, NSOL 

COMMuN/CnS  /SLLVEL»HOS  » VOS .COS,hPLS, VPOS,CPOS 

COMMjN/CNlC/XP0S»XSnL,VPSOL»VS0L  tCS  OL  »HU  S  t  VUS  »  Lo  S  »  MPUS  »  V  PU  S  .  C  P  Ui 
CUrtMON/Cnll/TL»TW,TT,SW 

C0HMGN/CK12,  XFAC2,KFAC  «  TAl^,f  hOV  ,f  ST  Al 

C 

G-9.3  I 
RHO-1000.0 
F8ICT»FSTAT 
FM1«RH0*C« ( HcS-TT/2.0 ) » T W« TT • T T /2. 0 
F«2-RH0«o«(TW2.O»»T'<»TT»0.5«(  HDS-TT»HDS» 
f M3-SK«C*Tl«C.5»C0S( ASIN( (HUS-TT)/TL) ) 
TX-(H0S-IT)/SIN(ASIN( ( HUS-TT  » / TL ) » 
fM^-XFAC2«P»-L»C*0.5*(HDS-TT)«TH«C.5*lX*«<: 
IF (HOS.LE .TT )  FN^-0.0 
rP«(0.5/TL)*lFN2*F13*FM'»-rNl) 

F  P?-(  FP-f         i0.b*JXU  ♦COS  (  ASIM  <HUS-I  T  »/TLn 
IUVSOL.GI.C.OFKICT-FNOV 
H5H-T  T»0.  b**!-.  PUS-IT  » 
CALL    SHAPt  <T  1  hE  ,HSrt.ASH,ISH,PEf.  ) 
FLOAk- ASM-T  T«  TH 
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CALL  SHAPC<TIHt»H«»US,ARUS«TUS,r£RI 
QSH-ARUS» (VPUS-VSJLI 
SUERV-OSH/ASM 

SHFRF-TAU«O.5«RH0»2.0»TL»Tl«SHE*lV*«2 
IFtHPUS.LT.TT }  Fr2-0.0 
FP2-XFAC«f Pi 

f SUR-FRIC7*<SW*C*C0S(ASINISO))-^P2l 
IF(FSUR.LT.O.O)  FSU^»0-0 

FS0L»RK0*C*T»»«0.5»<HPUS**2-«PDS«*2l»SK*C*SO-FSUR 
VPS0U-VSCL*DT*F5CL/SW 
hRITE<3«5  3)FS0L.FPZ.F  SUR 
53  FORHATI 10X,»F50L-» »Fl0.^t«   FP2-' «F 10. 4, •   FSUR=* vF 10.4,/ } 

IF<VPSOL.LE.O.Ot  VPSaL-0.0 
X!iOL-0.5»01«<  VPSOL»VSOLI*XSOL 
IFCVPSOL.Cl.O.O.ANO.VSOL.LE.0.0)  CCTO  5 
00  4  I-1»N 

IF<XNt  I  I.LE.XSOL.AMO.XMf  !<•  1 }  .  CT  .XSOLI  HSOL-I 

4  CONTINUE 

5  CONTINUE 
RETURN 
END 

c 
c 
c 
c 
c 
c 

SUBROUTINE  ASSIGN(N» 
C  THIS   SUBRCUTIhE    SETS   UP    THf    hEW   e*SE   CONOITIUNS   ALONO  THt 

C  PIPE    IN   FfFf-ARATION   f  Of   THE    ^EXT   Tint  SIcP. 

DIMENSION   OP  (61  )  ,0  rn(  30  )  ,T  l^(  3C  )  ,VP<  611  ,riP  (  61  )  tCP  (  fcl » 

DIMENSION  V(bl),H(6lt 

DIMENSION  VK  611  ,HR  Cijl  )  ,CS  <fci  )  ,VS(  »HS  (  61  l,CS<  61  )  ,SS  (61  ) 
DIMENSION  SP(fcl)»X'.(t>l}  ,XS  J6i  l,X^C61»  ,C(oi»  t 
COMM0^</CM3/NP^S»OIH♦^  IN 

C0HMGN/CM^/V,H,C»V-?,H^  ,Ck,XR,S^TVS,HS  tCS^XStSS^XN 
COMMON/tnt/CP ,VF,HP ,CP 
COMMON/CMfc/ iNSOLfNSHL l,NSOL 

C0MM0N/C«9/SCLVEL»  H3S  t  VDStCLS»HPl/S.  VPOStCP  OS 

C0MMCN/CM1C/XPCS,XSTL  «  VPSCL,VSCL  «OPSl)L»H.US»VUS»Cu  StHPUS,V?USt  CPUS 
00    1  I-^1,N«1 

IF  (VSOL.EC.C.C.ANO.I.  EQ.NSCL)    CCTO  2 

v(i)  =  vp(n 

H(I)«HP(I) 

c<  n-cp<  I » 

IF ( In  SOL .CT.O .AND. I .£ Q.NSCLl 1    COTO  6 
7  CONTINUE 

;F (VSCL.E C.C. C.ANO.NS OLl.CT.t )   COTC  2 
GOTO   1  *        '  .     •..  . 

2  HUS«MPUS 
VUS-VPUS 
CUS-CPUS 
HOS-HPOS 
VOS-VPOS 
COS-CPDS 
HPCNSOL )-HPUS 
VP  INSOD-VPUS 
CP(NSOL)-CPLS 
H(HSOL)»HPtS 
V»NSOL)-VPUS 
C<NSOLI«CPUS 

core  I 

6  HPUS-HP(JI 
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vpus-vpc I ) 

VPDS-VP<I > 
CPUS-CP( I ) 
CP0S-CP<1 J 
COTO  7 
1  CONTINUE 

SOLVEL-VSOL 
VSCL-VPSOL 

IFfVPSOL.CE.O.OI  GOTO  3 
GOTO  A 

3  VSCL-VPSOL 
HUS-HPOS 
»US-VPUS 
COS-CPUS 
HOS-HPDS 
VDS-VPOS 
COS-CPOS 

4  CONTINUE 
RETURN 
END 

SBENO 
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Appendix  2 


Deflvation  of  solid  boundary  equations 
used  in  Subroutine  SOLID 
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Derivation  of  solid  boundary  equations  employed  in  Subroutine  SOLID. 
Available  equations  at  moving  solid  boundary. 
1)    Leakage  flow  past  solid. 


Leakage 
flow,  Q 


Q  =  XK  {SE  -  SEq}2 

where  SE  is  the  specific  energy 
relative  to  moving  solid. 


9 

Specific  energy  h  +  V^^/Zg 


XK,  SEq  determined  by  measurement. 

2)    Flow  over  solid  in  terms  of  flow  area 


Vg  =  solid  velocity 


Relative  velocity  Vj-g  =       -  Vg 

Flow  area  A  =  f(h) 

Total  flow  =  V^A  =  V^fCh) 

Flow  past  solid  =  Vj-^A  =  Vj.gf(h). 

3)    Characteristic  equation,  C*"  in  both  sub  and  supercritical  flow  conditions. 

P 


At 
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C"^  links  R  at  t  to  P  at  t  +  At 
C"^  equation  Vp  =  X2  -  Zlhp 
XI  =  g/cR 

X2  =  Vr  +  ghR/cR  -  g(SR  -  So)At 

Note  charcterlstics  equations  in  terms  of  absolute  velocity. 
The  flow  Q  over  the  solid  is  given  by  the  empirical  relation,  ^ 

2 

Q  =  XK  (h  +  ^re  -  SE^)^ 
2g 

where  Vj.g  is  the  flow  velocity  relative  to  the  moving  solid. 

Q  =  V^gA  =  (Vabs  -  Vb)A  =  (X2  -  Xlh  -  Vb)A 

Also,  Q  =  XK  {SE  -  SEq}^ 

=  XK  {h  +  £^abs  -  ^B^^  -  SE^}  ^ 
2g 

Hence 

(X2  -  Xlh  -  Vg)A  =  XK{h  +  (X2  -  Xlh  -  Vg^^  -  SE^}^ 

2g 

Collecting  terms 

U  =  (X2  -  Vb)A/XK;  B  =  XI  A/XK; 

W  =  (Xl2)/2g;  Y  =  X22/2g  -  SE^  -  VgX2/g  +  V^/Zg; 
Z  =  Xl.vg/g  +  1.0  0  X2Xl/g 
U  -  Bh  =  {Zh  +  Y  +  Wh2}2 

=  z2h2  +  y2  +  w2h^  +  2YWh2  +  2ZYh  +  2ZWh3 
W2h^  +  2ZWh3  +  (z2  +  2YW)h2  +  (2ZY  +  B)h  +  (Y.lYl  -  U)  =  0 
Note:     Y.ABS   1y|  retains  sign  Y. 

This  quartic  is  solved  by  the  Newton  Raphson  method  as  B  and  U  involve  f(h) 
through  the  flow  area  terra. 
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APPENDIX  3 
Typical  Output  Program  TRANSCD 
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